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ABSTRACT 


The design of anti-icing protective devices for aircraft 
requires that they be checked for operational effectiveness in a testing 
facility. For one such facility, the icing wind tunnel, the method used 
to simulate the icing cloud is to inject a continuous spray of water 
droplets into a cold airstream. These simulated conditions were investi- 
gated in this thesis using a one-dimensional numerical model. 

The purpose of this study was to formulate a model of general 
applicability that determined the spray and air thermodynamics and 
dynamics in any test facility. The governing equations were programmed 
in FORTRAN IV language for solution on a digital computer. The temper- 
atures and velocities of the spray droplets and the air, averaged over a 
cross-sectional area perpendicular to the tunnel axis, were calculated. 
Two models -- a Single droplet model and a droplet ensemble model -- were 
developed. The equations of the two models treat a frictionless, 
adiabatic flow in which the total entropy of the system is conserved. 

The equations of the single droplet model apply to situations where the 
spray has a negligible influence on the air flow. The ensemble model is 
applicable for cases where the contributions of the spray (consisting of 
N identical droplets per unit volume) to the moist air have to be consid- 
ered. 

A limited parametric study was performed to determine the effects 
of moist air properties at the tunnel plenum (beginning of contraction) 


and of initial water spray conditions on the state of the flow at the 
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working section of the NRC* high speed tunnel. The results indi- 

cate that the downstream spray properties (at the tunnel working section 
section) )are essentially independent of the initial spray liquid 
content, injection temperature, and relative velocity (with respect 

to the air). The spray temperature is sensitive to the tunnel relative 
humidity and tends to approach the wet-bulb temperature of the air. 
Smaller droplets (10 microns) approach the wet-bulb temperature and the 
velocity of the air sooner than larger droplets (100 microns). When 
plenum conditions are adjusted to give air speeds at the working sec- 
tion of 30 m/s and 60 m/s, a spray of 25 micron droplets reaches the 
air temperature and velocity before arriving at the working section of 
the NRC tunnel. The spray remains substantially out of equilibrium with 


the air when the final air speed is near 130 nys. 


* National Research Council 
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CHAPTER 1 


INTRODUCTION 


1.1 An Introduction to the Icing Problem 

The formation of ice on aircraft surfaces occurs during flight 
through clouds of supercooled water droplets. Because these supercooled 
water droplets are in a state of metastable equilibrium, the disturbance 
created by a moving aircraft causes the impinging droplets to form ice on 
its nucleating surfaces. 

The ice accretion on aircraft surfaces always results ina 
degradation of its performance and reduces its operational safety. The 
result of ice accretion may be malfunctions of essential instruments, a 
reduction of aerodynamic lift, an increase in drag, an increase in weight, 
and an impairment of manoeuverability. Icing of propellers and rotor blades 
may also lead to a considerable reduction of forward thrust. In the case 
of turbine-type aircraft, ice released from the inlet duct of the engine 
may damage fan and compressor blades and may cause compressor flame-out 
(Vath, 1978). Because of these effects, aircraft are routinely equipped 
with protective systems for either removing accumulated ice or for 
continuously maintaining the exposed surfaces of the aircraft free from 
ice. These protective devices must be tested for operational effectiveness 
under conditions that closely simulate or duplicate the natural atmospheric 
icing conditions. A good simulation thus requires that the primary 
variables defining icing conditions be established and controlled in the 
testing facility during the testing process. The purpose of this thesis 


is to analyze the conditions that are simulated in an icing wind tunnel. 
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Here at the University of Alberta, the Department of Mech- 
anical Engineering has developed an icing wind tunnel to simulate nat- 
ural atmospheric icing conditions. An icing wind tunnel has also been 
developed at the National Research Council (NRC) in Ottawa. Chapter 4 
of this thesis describes the numerically sinulated icing conditions that 
result from a high speed air flow in the NRC tunnel. For a general des- 
cription of the NRC tunnel, the reader is referred to Section 1.3. The 


next section presents a survey of natural icing conditions. 


1.2 Natural Icing Conditions 

In general, ice formation results from a complex interaction 
of aerodynamic and meteorological factors. The aerodynamic factors include 
the aircraft speed, skin temperature, and the size and shape of the 
icing surfaces. The meteorological factors include pressure, air 
humidity, liquid water content and cloud ice, cloud droplet size and size 
distribution, and cloud temperature. The degree or severity of icing 
may be described as light, moderate, or severe (Vath, 1978), but these 
are relative terms. The difficulty with them, from a meteorologist's 
point of view, is that they are usually aircraft specific. Thus, light 
icing for a Boeing 747 may be severe for a helicopter, for example. 

There are three types of atmospheric ice formation, namely clear 
ice, rime ice, and mixed or intermediate ice. Clear ice generally forms 
with air temperatures in the range of 0°C to -10°C, where freezing is very 
slow. Here, larger supercooled droplets that collect on the aircraft 
surfaces partly freeze on impact, while the rest of the droplet runs back 


prior to freezing. The resulting ice formation is dense and adheres tightly 
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to the surface. The liquid water content of clouds producing clear icing 
is a large one (w1g/m>) because of the larger droplet size. Rime ice, 

on the other hand, forms at low temperatures (as low as -40°C) and usually 
forms on aircraft surfaces from flight in clouds with low liquid water 
contents. This ice form accumulates on the leading edges of aircraft 
wings and other surfaces and results from small droplets freezing on the 
airframe individually with little or no spreading. The physical appear- 
ance of rime ice, often described as friable, is milky white (due to the 
presence of air pockets) and its structure is porous and the surface 
rough. Rime ice tends to build forward from exposed surfaces and accum- 
ulates less quickly than clear ice due to the generally lower liquid 
water contents of clouds producing this type of icing. Mixed or inter- 
mediate ice, as its name suggests, has some characteristics of both types 
of ice. 

The three types of ice formation described above may be 
associated with the following major types of icing environments: 
stratiform or layer clouds, cumuliform clouds, and freezing rain or other 
low level icing. The stratiform clouds, with low to moderate liquid 
water contents, have a long horizontal extent and are associated with 
continuous icing. The cumuliform clouds, with moderate to high liquid 
water contents, have a vertical extent comparable to the horizontal 
extent, and are associated with intermittent icing. The freezing rain and 
freezing drizzle situations (i.e. moist layer(s) above the surface 
characterized by temperatures near freezing) are associated with severe 
continuous clear icing. 

The icing intensities associated with the different cloud 


types are described in Table 1.1. It is emphasized that various types 
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TABLE L. 


ICING POTENTIAL ASSOCIATED WITH THE VARIOUS 


CLOUD TYPES (AFTER MASON, 1957) 


Type of Cloud* Intensity of _ Probability that Water content 
Icing ice will accumulate (g/m) of cloud 
Tels CNS May be severe High 0.2 - 4.0 
SC,uAC 
Ac with As Light—Moderate 3503 Ole 0s 
As Light—Moderate Low O58 =2075 
Sc Light Low Oo =055 
= Tew Towering Cumulus 

Cb Cumulonimbus 

Ns Nimbostratus 

Sc Stratocumulus 

Ac Altocumulus 

As Altostratus 

Sst Stratus 


Cu Cumulus 
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and degrees of icing are possible in convective type clouds with clear 
icing often observed. In the case of Sc type clouds formed over the 

- sea and Sc formed by the spreading out of Cu, moderate to severe mixed 
or clear icing may occur. 

The fact that ice formation differs from aircraft type to 
aircraft type results in a rather subjective classification of the 
different degrees of icing. Despite this subjectivity, standard icing 
conditions have been formulated for the purpose of testing aircraft 
performance. The icing criteria are based on flights made with specially 
instrumented aircraft through weather conditions that have resulted in 
icing. Since these standard icing conditions are defined in terms of 
the typical speeds, rates of climb and performance of a particular 
aviation era, the standard icing conditions can become outdated as air- 
craft performance improves. It is emphasized that the specifications 
used in testing an aircraft's ability to withstand icing pertain specif- 
ically to icing on external aerodynamic surfaces and not to engine 


performance icing conditions. 


1.3 The NRC High Speed Icing Tunnel 

A vertical cross-section of the NRC high speed icing wind 
tunnel is shown in Figure 1.1. The portion of the tunnel where icing 
conditions are simulated consists of 1) a rapidly converging bell- 
mouth contraction followed by 2) a more gentle linear contraction and 
3) a working section of constant dimensions. (See upper left hand por- 
tion of Figure 1.1.) An analytic expression for the variation of the 
cross-sectional area along the tunnel axis appears in Appendix H. This 


cross-section is square throughout the tunnel. 
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Briefly, the method used at NRC to produce an icing cloud 
is to inject a continuous spray of water droplets into a cold airstream. 
The droplet spray is produced by four pneumatic atomizing nozzles 
located 0.20 m downstream from the plenum (beginning of contraction). 
The water supply to the spray nozzles is usually at room temperature 
while the compressed air supply needed to produce the spray is heated 
to 30°or 40°C. Since both the water flow rate and the atomizing air 
pressure are roughly controllable, the liquid water content of the spray 
and median volume radius of the droplets are roughly known. On the other 
hand, the spray temperature, injection velocity, and the precise location 
of the droplet formation are unknown. 

For most experimental and modeling situations, it is assumed 
that the droplet attains the temperature and velocity of the airflow 
by the time that the spray reaches the working or test section of the 
tunnel. This assumption is generally valid for slower tunnel speeds, 
For the higher tunnel speeds (greater than 90 m/s in the working section) 
and large accelerations, there may be insufficient time for the spray to 
reach thermal and dynamical equilibrium with the airstream, especially 
if the airstream itself is cooling by adiabatic expansion and acceler- 


ating. 


1.4 Purpose of the Study 

As mentioned earlier, the purpose of this study is to analyze 
the conditions that are simulated in an icing tunnel. The problem to be 
considered is to determine the temperature and velocity of the spray and 
moist air in the working section of the tunnel. For many testing situa- 


tions, measurements of these properties are either not feasible or are 
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subject to significant error, It is also desirable to determine the 
initial conditions that are required in order to produce given icing 
conditions at the working section. The approach used in this study is to 
numerically determine these flow properties from a model. The numerical 
models that were developed to simulate these conditions are discussed in 
Chapters 2 and 3. Results describing the numerical simulation in the 


NRC high speed tunnel are discussed in Chapter 4. 
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CHAPTER 2 


THE SINGLE DROPLET MODEL 


2.1 Outline 

The liquid water contents found in a cloud environment are 
simulated in an icing tunnel by injecting water droplets, known as a 
spray, from atomizing nozzles into a moist air flow. The problem of 
modeling the thermodynamics of the injected spray is to describe the two 
phase flow of a water suspension in air. The medium considered consists 
of moist air and fine liquid droplets that are, for practical purposes, 
accelerated by the viscous drag of the accelerating air. The droplets 
are assumed to be spherical, rigid, and incompressible and their specific 
volume is assumed negligible compared to the gas phase. In addition, the 
effects of droplet coagulation are neglected. 

A description of the dynamics of the two phase flow uses the 
fundamental conservation laws for mass, energy, and momentum. The 
balance equations are written for air and for the droplet, respectively. 
These equations account for mass, momentum, and energy exchanges between 
the air and the droplets. 

The two sets of equations, with their inherent assumptions, 
constitute the models of the present study. The equations are solved 
numerically. One set of equations describes the properties of the system 
where a single droplet is injected into the moist air flow. The second 
set describes the properties of the system where an ensemble of N identical 
droplets per unit volume is injected into the air flow. In the discus- 
sions that follow, the first model is referred to as the "single droplet 


model" while the second model is referred to as the "ensemble model". 
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Both models treat the flow in the tunnel as a one-dimensional steady 
flow. 

The two-dimensional nature of the tunnel flow is ignored in 
both models. Instead, general one-dimensional equations are obtained 
from the use of spatial averaging techniques. The variables used in 
this study are considered to be averaged over a cross-sectional area 
that is normal to the tunnel axis and are functions only of the distance 
x along the tunnel axis from the plenum (beginning of contraction). 

A number of other averaging techniques may be employed in engineering 
applications when treating one-dimensional averaged variables. These 
include (Solbrig and Hughes, 1978, p. 310) volume averaging of variables, 
spatial averaging of the equations of motion, and statistical methods of 
averaging. 

The advantage of using one-dimensional averaged variables is 
that the resulting flow equations are simplified and become very general. 
When using averaged variables, however, it is important that the infor- 
mation obtained from solving the flow equations be carefully interpreted. 

Appendix 1 contains a summary of the variables used in the text. 
Variables are also defined in the text as they are introduced. Variables 
describing properties of the dry air, vapor, and the droplets are denoted 
by the subscripts d, v, and c (condensed phase), respectively, while 
variables without subscripts refer to properties of the moist air (dry 


air plus water vapor). 


2.2.4 Dry, AimeFilow: inva Wind Tunnel: 
Many wind tunnels are designed with an objective of obtaining 
an air flow in the working section that has a uniform, rectilinear 


velocity profile. For a properly designed tunnel, the velocity over a 
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cross-section of the working section (outside of a thin boundary layer), 
described as a flow area, may be treated as constant. In the contraction 
section, the camponent of velocity perpendicular to the tunnel axis is 
treated as negligible if the curvature of the locus of the center points 
of the flow areas is small and if the tunnel boundary layer is thin. In 
such cases, the flow area is approximately the cross-sectional area and 
the averaged velocity is approximately the axial velocity. 

The following analysis of one-dimensional flow in a wind tunnel 
describes a positive distance x as the distance along the tunnel axis 
fran an origin at the plenum. The portion of the tunnel considered is 
the contraction and working sections. 

The relationship between pressure and velocity at different 
points of a fluid field can, in general, be obtained from an analysis 
of the dynamics of a particle of the fluid. The relation is obtained 
from the simultaneous solution of the energy, continuity, and momentum 
equations, and from the equation of state. 

The equation of state adopted for the dry air is the Perfect 


Gas Law, written as 


Cae foe Demta} 
ea hy 


where 

P is pressure, 

0 is density, 

T is temperature, and 

R, is the specific gas constant for dry air. 


Corrections for moist air, to avoid a variable gas constant, are made 


later. 
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The momentum equations describing the motion are the Navier- 
Stokes equations. Their complete solution is very complex unless a number 
of simplifying assumptions can be applied to the flow. The equations of 
motion reduce to the compressible Bernoulli equation if it is assumed 
that the flow is: 
1) steady (does not vary with time locally) and, 
2) inviscid (viscous or frictional forces are negligible), 
If the flow is treated as a one-dimensional flow, the analysis is 
simplified further. The Bernoulli equation for a one-dimensional flow 
can be written as: 


P 
U2 yu. 2 
2 i 

Se 


git 2a) =" 0 (2 22) 


omy 2 
where U refers to velocity, z is the vertical distance from the tunnel 
axis, and g is the acceleration of gravity. The subscripts 1 and 2 
refer to values of the variables at the positions x) and X51 respectively. 
For the tunnel with velocities greater than 10 m/s and vertical scale of a 
few meters, scale analysis shows that the gravitational term in the above 
expression is neglected relative to the other terms. 

When the walls of the tunnel are assumed to be insulated and 
all other heat sources or sinks are negligible, the steady-state thermo- 
dynamic equation describes an isentropic flow of a perfect gas, The 


adiabatic expansion and compression of the air is described by the rela- 


tion; 


= oe (2)..2 3) 


and y = ee =) 4¢for dry*air. 


pee eigen Birds etl 
if ermal (is) sts aCe pret 
Peri aeasir “i Ny, iad eatogaiiag ite se os 


So sola ot pas ot ta 
homewar si nest ee is iadae 


. 


: Soe vibeped rr ibe 


ASI Pal Drei ets 96 pups 


oh) tunaet tare eat’ ta ae i <i 


be ‘ones oe ae Pcie) PE) Ls i a a ad ‘i 
ay. i 


aoe Tan ee meoriete we be ‘EAS ek s dab oy athe 


Site t aaainpehls OO. 2 vie Ne socom silt 
Yea taper: a gh AEST See MNT, Be ian \k 


b-99 @LADe lepisagy bi: 
4 | ) ad 

VER eS! cu eS adh sey dail co: canker ofp egies aes 

} ater 


: aamtis. tel Sic ei et wher ® | 
ry Ce ce ae ae os "Bova "SE Feayenid List ‘inition a 


“ie 


_ or wis 7* wie ® e¢*y - ar Riel oar pe vai * aries eae 


7 nen 


- “a | gate sath Aree i iy Sir c}- thes! peta i aang 


Allens tity col teks: ai te “igh, me Tas te bcd iti eal 


13 


Finally, if mass is conserved, the steady-state mass continuity 
equation is easily integrated for the one-dimensional flow. The equation 


describing mass continuity of the dry air flux is, 


le) UnA, = p, UA, (2e2Fy) 


where A is the cross-sectional area. 

For adiabatic flow where mass is conserved, the momentum and 
energy equations can be combined to give an expression for the velocity. 
The relations describing the state of the moving air are as follows: 

1. Bernoulli equation (2.2.2), 

2. Dry adiabatic expansion (2.2.3), 

3. Mass continuity equation (2.2.4). 

If (2.2.3) is substituted into (2.2.2), neglecting the gravity term, the 


resulting equation may be integrated to give, 
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2 y= 2 2 Ged l 


Equations 2.2.3 and 2.2.4 are combined to give the expression, 


v=] 
P P UA 
Page Biel 1 ’ (25266) 
Py | UA, 

Tf (2.2.6).is| substituted:into (2.2.5), an-equation.is obtained jthat 
relates the fluid velocity and the cross-sectional area at a point along 
the tunnel axis to known variables at another point. The equation that 


results, 
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can then be solved numerically to yield the velocity, Us. for any cross- 
sectional area in the x-direction. 
Using the one-dimensional velocity, the density is determined 


from the continuity equation using the relation, 


The remaining equations, the adiabatic expansion relation and the 

Perfect Gas Law, can be applied to obtain the pressure and the temperature, 
respectively. In addition, the dimensionless Mach number, M, can be 
calculated using the above quantities. The Mach number for adiabatic 


flow of the dry air (treated as a perfect gas) is given by, 


Mriesies eu eet oo 2.8) 


Equations 2.2.2, 2.2.3, and 2.2.4 can also be rewritten in terms of the 
Mach number. The derivations are given in Appendix 2. 

A numerical solution to (2.2.7) can be obtained using the 
observation that us/2 is small compared to the second term on the left 
hand side for velocities of interest. If (2.2.7) is rewritten as 
| y2 


U2 ey a, 
pelrdmitles | its ees ede ee a (2.2509) 
co Py A 2 2 2 Py yal 


it is found that the resulting expression is particularly suited to 
yield a numerical solution for the velocity by an iterative technique. 
A first guess of the unknown velocity, Us, on the right hand side is 
the velocity obtained under the assumption of incompressible flow. 


The expression for the first guess velocity, denoted Us ie is 
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An improved estimate of the root is written as, 


Can Ney, 
Gen! 2 Qane ic 2P)  [UiA) ) 
U = = as Mee x rs 
(U, Uo) + iz a / a es | 


(222.10) 


(n + 1) 


where U refers to the (n + 1)th estimate of the one-dimensional 


velocity, U,. 


2.3 Moist Air Flow 

When modeling the moist air flow in an icing tunnel, the 
equations of motion must account for the presence of water vapor. In 
practice, the water vapor contribution to the moist air density is 
relatively small and, depending on the accuracy required, some calcula- 
tions may neglect this component. In the following discussion, the vapor 
contribution to the moist air density is treated as constant since it is 
assumed that evaporation or condensation does not occur. As a result, 
the inclusion of water vapor will alter only the definitions of the 
variables characterizing the flow. 

A number of variables are defined in the meteorological liter- 
ature to describe the water vapor content of moist air. One such moisture 
variable, the mixing ratio, r, is defined as the ratio of the density of 
water vapor to the density of dry air; that is, r= p/P . Using the 
Perfect Gas Law for each component, this ratio is also written as, r = Dp 
where ¢ = RJR, = 0.62198 and R, is the specific gas constant for water 
vapor. P. and e are the dry air and water vapor partial pressures, 


d 
respectively. 
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A humidity variable, virtual temperature, incorporates a 
moisture correction into the temperature. The virtual temperature is 
_ defined as the temperature at which dry air would have the same density as 
a sample of moist air, assuming both have the same pressure. Its use 
allows the equation of state for dry air to be written for the moist 
air without resorting to a humidity-dependent gas constant. The general 


formula for the wirtual temperature, T , is 


(2e3e)) 


Tl’ =T \ 1 eT 
r+] 


Using the definition of virtual temperature, the equation of state for 


moist air is, 


P = = - 
zs +e p RY u 


where p = Pg + Py When reference is made to the moist air temperature 
(which for a dry air and vapor system in thermal equilibrium is also 
the temperature of the dry air and of the vapor), it will be denoted as T. 
As mentioned, only the definitions of the parameters and vari- 
ables used in the equations of motion for the dry air flow will require 
modification when water vapor is a part of the air flow. For example, 
the equation of state for moist air specifies that the pressure and 
density refer to their moist air values and that the temperature used 
is the virtual temperature. The Bernoulli equation is likewise unchanged 
in form when the moist air variables are used and is given by (2.2.2), 
Since it is assumed that the mixing ratio remains constant for this flow, 
the steady-state continuity equation for the moist air implies continuity 
of the dry air flux, Equation 2.2.4. In addition, Poisson's equation 


(2.2.3) for an adiabatic expansion of an ideal gas from the state Pi rPy 
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to the state P5105 applies for the moist air flow if the exponent, 
y= Ss , is determined as a function of the mixing ratio of the moist 
air. Expressions for the heat capacities of moist air, and Cc, , are 


derived in Iribarne and Godson (1973, p. 70) and can be written as, 


Cc Cc 
r pv f VV 
Cg ec 1+—~(—--1)/ , cc =c 1+—(—- 1) 
p pd fal Sari | Vv vd | rr Cg 
(23. 2,) 


Here, Co and Cra are the specific heat capacities at constant pressure 
for water vapor and dry air, respectively and one. and Cod refer to 
specific heat capacities at constant volume for water vapor and dry air, 
respectively. It is seen from (2.3.2) that the specific heat, 
capacities remain constant for a constant mixing ratio, as does the 
exponent y . 

Since the equations describing the moist air flow are identical 
in form to those describing the dry air flow when the mixing ratio is 
constant, the velocity equation derived in the previous section applies 
to the moist air flow. The one-dimensional velocity is then solved 
using the iteration equation (2.2.10). Following the analysis of the 
previous section, the moist air density is determined from the continuity 
equation. The moist air pressure, virtual temperature, and the thermo- 
dynamic temperature, T, are then found using the adiabatic expansion 


equation, the Perfect Gas Law, and the definition of virtual temperature, 


respectively. 


2.4 Droplet Motion 


The motion of a single droplet in air is governed by Newton's 
Second Law of Motion, where the forces considered are the gravity force 


and the drag of the air on the droplet. The gravity force may be neglected 
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relative to the drag force for the droplet sizes (of the order of 
10-100 microns) and accelerations considered in this study. 

The problem of determining the droplet motion in a viscous 
medium is a camplicated one since the drag force is a function of the 
droplet size, relative velocity of the droplet with respect to the air, 
the viscosity of the air, and of time. Three of these variables can be 
combined to form a ratio, the Reynolds number, that characterizes the 


viscous flow. The ratio is given by, 


2r i Hira UN U | 


= = : Zea | 
fe. = lay UN : ( ) 


| where the variable, Loe refers to the radius of the droplet, u is the 
dynamic viscosity of the air, and v = p/p is the kinematic viscosity of 
the air. It is assumed that the droplet remains spherical, 

A related problem, the determination of the steady-state drag 
on a sphere, has been the subject of a considerable number of theoretical 
and empirical studies. A convenient summary of these investigations, as 
well as a discussion of a numerical study to extend the applicable range 
of the past investigations, is presented in a paper by IeClair et al. (1970). 

For droplets approximated as rigid spheres, the numerical study 
of IeClair et al.(1970) is able to theoretically predict accurate values 
of the drag over a Reynolds number range from 0.01 to 400. Their theor- 
etically determined solutions agree excellently with experimental drag 
determinations by Beard and Pruppacher (1969). The study by Beard and 
Pruppacher measures the drag on water drops falling at terminal velocity 
in air for the Reynolds number range 0.2 to 200 and extrapolates these 
results to a Reynolds number of 400, numerically. The agreement of the 


theoretical and numerical results is given as evidence that the droplets 
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can be modeled as rigid spheres for Reynolds numbers in this range, The 
empirical drag relation described by Beard and Pruppacher's results is, 
F 


Csasies 


= 1 + ake (22422) 


where 


&@= 0.102, 8 = 0.995 for 0.01 = Re = 1,5, 


Il 
1A 


202, 6 Ges02e fOr, Leo Re < 20, 


Re < 400 


LA 


a= 0.189, 8 


II 


O2052 ,for 20 


Here, F is the steady-state drag on the spherical droplet and Sioa pel is 


the Stokes drag given by, 


iotokes * 6rur (U . Uo) 


The steady-state drag given by Equation 2.4.2 is not adequate 
for many situations where the droplet undergoes large accelerations. For 
example, preliminary calculations for droplet motion through the contraction 
section of the high speed icing tunnel at NRC in Ottawa (to be discussed 
later) indicate that additional terms should be included in the equations 
of motion to account for the time dependency. An analytical solution to 
the non-steady state problem, for Reynolds numbers less than 0.01, has 
been obtained by Pearcy and Hill (1956). The drag force predicted by 


solving the Navier-Stokes equations is 
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The first term is the steady-state viscous force arising from 
uniform motion (i.e. Stokes drag). The second force term is associated 
with the momentum transferred to the surrounding medium by changes in the 
velocity of the droplet. The third term, known as the history term, 
represents resistance due to the expenditure of energy needed to set the 
fluid in motion. The integrand of the history term indicates that the 


1/2) 


effect of past accelerations decreases with time (proportional to t~ 
The Stokes drag term in Equation 2.4.3 is corrected using the 


steady-state results of Beard and Pruppacher (Equation 2.4.2). This 


leads to | 
du 3u U-U lad 
Cush 3 Cc B a ene A 
a Ems 19 ote Seen (ist GRE 4 (U-U_) 
drag Mo dt re { ae ( Ete Cc 
Cc 
Set di U-U) 
Be 838 ‘Pp rs ais Cae \ (2.4.4) 
G mL dt ie 


In most of the literature on droplet motion, the steady-state 
drag forces are used to predict the motion of water droplets in a viscous 
medium. Justification for neglecting the history term is generally not 
offered. An empirical study by Sartor and Abbott (1975) concludes that 
steady-state drag forces can be used to describe the motion of an accel- 
erating droplet that is falling in a straight path and is characterized 
by a Reynolds number less than 5. An investigation by P. I. Joe (1975) 
highly recommends that the history term be included in the drag formu- 
lation of accelerating droplets when the motion is along curvilinear 
paths and for large Reynolds numbers. The study shows that, for inter- 
mediate or large Reynolds numbers (Re > 10), errors of at least 25 per- 


cent can result when calculating droplet velocities. 
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If the full history term is included, determination of the 
time dependent motion is not trivial, since the solution of the differ- 
ential equation for the droplet motion (Equation 2.4.4) first requires 
integration of the history term. In order to simplify the calculations, 
the history term is omitted from the drag equation that is solved 
numerically in both models. It is realized that the error due to its 
neglect may be significant. 

Figures 2.la and 2.1b show the results of a numerical exper- 
iment performed to estimate the relative importance of the history term 
for high speed flow in the NRC tunnel. The experiment calculates the 
ratio of the force given by the history term to the steady-state hori- 
zontal drag force as a function of the distance, x, the droplet travel 
time, and the droplet initial conditions. A first approximation to the 
history term is obtained using velocities that were given by the equa- 
tions of motion from which the history term had been omitted. The drop- 
let radius used in this experiment is 10 microns. The air speed at the 
sprayers is 51 m/s and the corresponding air temperature LS eG; 

The Figures indicate that the relative importance of the 
history term is strongly dependent upon the droplet initial conditions 
(that is, the droplet history). It is seen that the drag force is always 
of greater magnitude than the history term. For the case where the 
droplet injection velocity is zero (and the initial relative velocities 
are large), the Figure indicates that the history term is relatively 
unimportant for a large fraction of the travel time. It is noted, how- 
ever, that the history term is of significance for droplet motion near 
the sprayers. For the case where the droplet is injected with the air 


speed, the initial drag forces are extremely small and the history term 
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FIGURE 2.la (TOP) and FIGURE 2.1b (BOTTOM) : Importance of the 
history term relative to the drag force as a function of 
injection velocity. The arrow in Figure 2.1b indicates 
the beginning of the working section. 
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approaches a maximum of relative importance soon after injection from 
the sprayers. When the droplet is near the sprayers, the history term 
should not be neglected if the initial injection relative velocities 


are small. 


2.5 Droplet Growth and Thermodynamics 

When a pure water droplet is placed in an environment of humid 
(or dry) air, the droplet will grow or evaporate by the net diffusion 
of water vapor molecules towards or away from its surface, Associated 
with this diffusional growth is a substantial flow of heat to or from 
the droplet surface as a result of the release or absorption of latent 
heat. 

Any attempt to solve the problem of diffusional growth by 
condensation must consider the droplet thermodynamic and growth equations 
Simultaneously. The simplest case, which is treated in this section, 
describes thé processes that control the growth of a droplet that is 
stationary with respect to the surrounding air. In the following section, 
the growth of a moving droplet is considered. 

The assumptions considered in the single droplet and the droplet 
ensemble models are similar to those used by many authors in modeling 
diffusional growth. The simplified growth and thermodynamic equations 
that result have been derived by Pruppacher and Klett (1977), for example, 
In short, the principal assumption is radially symmetrical diffusion of 
vapor to or from a motionless droplet. Additional assumptions that are 
used include the eelieiic: 

1. The gas that is immediately in contact with the droplet surface is 


in equilibrium with the droplet. Hence, the temperature and vapor 
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density fields immediately over the droplet surface are determined 
by the properties of the droplet. 
2. The steady-state equations are used, 
3. The droplet can be modeled as a rigid sphere of uniform temperature. 
A further discussion of these assumptions and their validity for the 
present study is given in Section 2.7. 

The forms of energy accounted for in this model include latent 
heat, the internal energy of the droplet, radiation, heat conducted to 
or from the droplet surface, and convection. The last term refers to 
the effect of droplet motion on the heat transfer and is discussed in the 
next section. 

The rate of droplet growth in air is governed by the rate of 
transfer of water vapor between the droplet surface and the environmental 
air. Using the above assumptions, Pruppacher and Klett (1978, p, 414) 
formulate an equation for the mass rate of change of the droplet, dm_/dt, 


due to mass diffusion. The steady-state result is, 


Me = Amr 2p oS = dnd, (0) - 0 ) (25a) 


dt c 


where & is the diffusivity of water vapor in air and p,, and p,, refer to 
Cc 


vapor density fields at the droplet surface and in the distant (undisturbed) 
environment, respectively. It is assumed that the droplets are suffic- 


iently far apart so that it is meaningful to refer to a distant undisturbed 


environment. 


The latent heat, q, released by condensation at the droplet 


surface is given by 


dq : 
dt BUR AALS onl (2.5.2) 
latent c 


heat 
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where the specific latent heat of vaporization is denoted as hoe If the 
Perfect Gas Law is applied, Equation 2.5.2 can be written as, 


dq : eel * 


dt RT (2'..5,.13)) 


latent 
heat 


Following the rotation of the previous sections, Qe. is the partial vapor 
pressure at the droplet surface and e is the partial vapor pressure of 
the undisturbed environment. 

An expression for the process of thermal diffusion is formulated 
by Pruppacher and Klett (1978, p. 418) using mathematical arguments that 
are similar to those required to derive (2.5.1). The result describing 
the rate of heat loss by conduction is, 
oe her kisah) (2.5.4) 

conduction on dn 
where k is the coefficient of thermal conductivity of moist air, 

The heat loss by long wave radiation is obtained from the 
Stefan-Boltzmann Law (see, for example, Kern, 1950, p. 74) and is given by 
dq = hor? oE(T + - T*) Wales) 

: radiation - z 
where o is the Stefan-Boltzmann constant and E is the emissivity of the 
water droplet. The effect of short wave radiation on the droplet thermo- 
dynamics is ignored. 

The rate of change of the internal energy is written as, 


dT dr 
dq = e mrp Cat Ane 0 Tore ae (2.526) 
W 
«| internal ee ae at ae 


where Cc. is the specific heat capacity of liquid water. 
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Combining these terms, the energy balance equation for a stat- 
ionary droplet growing in a motionless atmosphere is found fram a balance 
of the internal energy, the latent heat, and the heat radiated and cond- 


ucted fram the droplet surface. The resulting heat balance equation is, 


4 dT hor Q AY 
ei CMe) = 
3 Co WrW a Oe Ra 


Eee pac 
+ | + gate Lhe) 


7h ale weed) Ty) (5 a7s 


where the second term on the right hand side of (2.5.6) is neglected 
relative to the latent heat term, ae 

The present study assumes that the pure water droplet exists in 
a supercooled state for droplet temperatures below 0°C. Experiments con- 
firm, however, that there is a limit to the degree of this supercooling. 
Factors that affect the maximum supercooling include droplet size, the 
rate of cooling of the droplet, and the concentration of impurities 
(Willbanks, 1973, p.14). When droplet freeze-out (i.e. freezing of the 
droplet) occurs, the thermodynamic and growth equations must account for 
the changes that freeze-out produces in the specific heats, the energy 
flows, and in the liquid water contents (See Section 3.1). Ina real 
test facility, droplet freeze-out may occur due to a contamination of 
the spray water with particles that cause ice nucleation. For the pres- 
ent study, the problem of droplet freeze-out is not considered. Given 
further continuation of this work, it is recommended that droplet freeze- 
out be included. 

Before the droplet growth and thermodynamic equations can be 
integrated, it is necessary to write an expression for the equilibrium 


vapor pressure, @., at the droplet surface. (It is noted that the droplet 
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would not grow if the vapor pressure actually attained its equilibrium 
value.) An appropriate expression for vapor pressure is Kelvin's 
equation (see Neiburger and Chien, 1960). This equation accounts for 
the effect of the droplet curvature on the vapor pressure of a pure 


water droplet and is written as, 


ATS 
ep = Suh) exp Gere Rate (2.528) 


The quantity o” refers to the surface tension of water in air and 


W/A 
Co5¢ (1) is the saturation vapor pressure with respect to a plane water 
surface. Both of these quantities may be calculated as functions of 


temperature. The quantity, e (T), 1s obtained from the Clausius- 


sat 


Clapeyron equation, which is discussed in Appendix D. 


2.6 Ventilation Effects 

If the air and the droplet are moving with respect to each 
other, the transfer of water vapor and heat between the droplet surface 
and the air is significantly enhanced by the motion of the air around the 
droplet. This ventilation effect results because parcels of air which 
are approaching the droplet temperature and vapor pressure by contact 
with the droplet surface are being removed continuously. 

The customary approach to this problem is to introduce dimen- 
sionless quantities that characterize the forced convection into the 
expressions for heat and mass transfer. One quantity that is commonly 
used in the chemical engineering literature is a mean Sherwood number, 
Sh. In the meteorological literature, the effect of ventilation is gen- 


erally expressed by means of a mean ventilation coefficient, £, defined by 
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dm 4 dm 
ae i ae} (2.6.1) 


where dm_/at is the evaporation rate of a moving droplet and (am_/at) 0 
is the evaporation rate of a stationary droplet. Since the mass flux 
will not be uniform over the ventilated droplet surface, a mean parameter 
f£ is used to characterize mass transfer from the droplet as a whole. The 


dimensionless quantities, £f and Sh, are related by the equation, 


Sh 


Theoretical studies have had limitéd success in describing the 
effect of forced convection on the heat and mass transport from a spherical 
body. As a result, the relevant equations (Navier-Stokes, growth and 
energy equations) are usually solved by finite-difference techniques. 

A numerical solution to the problem is presented in a paper by Woo and 
Hamielec (1971), for example. These results are in excellent agreement 
with wind tunnel measurements of droplet evaporation by Beard and Pruppacher 
(1971). The values of f that were obtained by Beard and Pruppacher are 


closely approximated by the following empirical expression: 


ee 


ae = J 400 8 of oRkiec REZ), 
for Scl73 Rel/2 < 1.4 
and, 
fF = a = 0.78 + 0.308 (Sc!/3 Re!/2), 
for Scl/3 Re!/2 > 1.4 (2.6.2) 
The dimensionless ratios used in (2.6.2) are the Schmidt number, Sc = v&y, 


and the droplet Reynolds number, Re. The last relation is valid for 
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scl/3p.l/2 < 16, which roughly corresponds to flow with Re < 320. 


The conventions that exist in the literature for describing the 
. Sherwood number are often ambiguous since many of these definitions are 
based on experimental or theoretical convenience. The Sherwood number 
used in this study is defined by Beard and Pruppacher on the basis of 
experimental convenience and is given by the relation, 


dm. i ine M¢ (e -e) Sh (2.6.3) 
tao i-ec/P a 2 


The above expression introduces the quantities, 


Ler F e.* e 
ractreumien ince septs, #5. eee) 


Which are respectively the mean or "effective" values over a vapor trans- 
fer path of the temperature, vapor pressure and diffusivity of water 
vapor. 

The effect of ventilation on the transfer of heat from the 
droplet surface by conduction must be considered. By analogy with (2.6.2) 
a dimensionless number characterizing the heat transfer , the mean 
Nusselt number, Nu, is introduced. The rate of heat loss due to conduc- 
tion and convection is given by 


Nu 


dq = nr k(T - T) (2.6.4) 


iE e 
conduction 
and convection 


The expression for the mean Nusselt number is similar to the 
expression for the mean Sherwood number. The mean Nusselt number is 
written in terms of the Reynolds number and the Prandtl number, The 


Prandtl number, Pr = v/«, is analagous to the Schmidt number and is 
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obtained by replacing the diffusivity of water vapor in the Schmidt 
number by the thermal diffusivity of air, K = Kee): The mean Nusselt 

Pons 
—+ = 1.00 + 0.108 (Pr¥3 Rel/2) , 
for Pri/3 Rel/2 eet 


and, 
ms = 0.78 + 0.308 (Pri/3 Rel/2), 


oem pa 2 eh G70.) 


The last relation is valid for Pe < 16 or roughly for flow with 
Re < 320. 

iE ene rerLeets Of venti lativon described by (2,6,3) and. (2.6.4) 
are introduced into Equation 2.5.7, a thermodynamic equation is obtained 


for the moving droplet. The relation is, 


dT hor 2 off el=< a 
aes ee Cae Toh CVs heen cut meh _7 ) Nu 
3 Be Ca ae R E | ome Aor K(T 1) 5 

oe AnoEré (To - o) (2.6.5) 


The growth equation for the moving droplet is given by (2.6.3). Using 
(2.4.4), (2.6.3), and (2.6.6), the system of equations describing droplet 


motion in an accelerating air flow can then be solved numerically. 


2.7 Discussion of Assumptions 

Tf the results of any model are to be useful, the assumptions 
that are applied must be realistic. With this objective in mind, this 
section evaluates the validity of some of the assumptions that are applied 
to this model. The same assumptions that describe droplet behavior are 


used in the ensemble model. 
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It was mentioned in Section 2,4 that significant errors may 
result when steady-state relations are used to predict the motion of a 
droplet that accelerates rapidly in a high speed air flow. For such 
cases, a number of non-steady state processes (that include the time 
dependent history term) should be considered. For example, observa- 
tions indicate that the drag of the air on large water drops, character- 
ized by a Re 2 300, progressively increases with Reynolds number above 
that for rigid spheres. This increase in drag is attributed to three 
processes -- oscillation of the drops, generation of internal circulations, 
and distortion of the drop shape from spherical. These processes, which 
are important for large drops, may also be important for smaller droplets 
when their accelerations (and relative velocities) are large. 

The assumption that water droplets can be modeled as rigid 
spheres is justified for steady-state motion with Re = 200 (LeClair et al., 
1969). This corresponds roughly to droplets with radii less than 
450 microns. For non-steady state motion (large droplet accelerations) , 
it is expected that greater surface stresses will result in a flattening 
of the droplet. As a result, the properties predicted by the drag, 
growth, and thermodynamic equations, which model the droplet as a rigid 
sphere, become inaccurate. 

The existence of an internal circulation inside falling drops 
is supported by many experiments. For a complete description of this 
effect and the effect of oscillations on drop. motion, the reader is 
referred to Pruppacher and Klett (1978, p. 305-328). In short, past 
investigations indicate that the contribution of internal circulations 
to the increase in the drag is negligible compared to the contribution 


from distortion in the drop shape. The importance of vibration to the 
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drag depends upon its effect on the drop shape. 

Detailed discussions of the assumptions needed to derive the 
steady-state growth and thermodynamic equations can be found in several 
studies. Some conclusions based on these studies are considered below. 
1. The time dependent solution to the problem of water vapor diffusion 

to a sphere (Borovikov et al, 1961) shows that the steady-state 
diffusional growth relation (2.5.2) can be used when the droplets 
are smaller than 100 microns in radius. Pruppacher and Klett (1978, 
p. 414) calculate that the steady-state description is not valid 
when significant changes in the ambient vapor field occur during 


droplet travel times less than a critical time, t 


= / (1) « 


crit 


The time, t , 1S an estimate of the time that is required to 


Erith 
establish steady-state diffusion if the droplet is growing in an 
undisturbed environment. For T = 0°C and P = 1 atmosphere, the 


calculations indicate that t < 1.7 microseconds for droplets with 


crite 
a radius smaller than 100 microns. Since the thermal diffusion 
problem is analagous to the mass diffusion problem and k#D, 
Pruppacher and Klett estimate that the steady-state thermal relations 
may also be used if significant changes in the ambient air temper- 
ature do not occur during droplet travel times less than torit: 
For the present study, changes in the ambient air proper- 
ties occur on a time scale of milliseconds so that steady-state 
relations may be used. 
2. It is customary to assume that the droplet has an isothermal surface 
because of the large thermal conductivity of water and because it is 


expected that small differences will exist in the heat flux around 


the drop surface (Ibid., p. 444). Contrary evidence from studies 
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showing that the ventilation effect has a strong local variation 
over the droplet surface suggest, however, that there must be a non- 
zero temperature gradient at the droplet surface, although the gradient 
will be somewhat degraded by the droplet internal circulation. 
For lack of information on this point, it is usually assumed that the 
droplet has an isothermal surface. 
The assumption that the moist air is a continuous field right up to 
the droplet surface is unrealistic for droplets with radii comparable 
to a mean free path of vapor molecules in air. The molecular exchange 
processes that occur at the liquid-vapor surface can be incorporated 
into a modified growth theory (Fukuta and Walter, 1970) using results 
from kinetic theory. Since the modified equations become complicated 
to parameterize and to solve, these molecular effects are not consid- 
ered in the present study. It is emphasized though, that errors may 
be introduced when the molecular processes are neglected. A study 
by Fitzgerald (1972, p. 132), for example, calculates that the 
instantaneous growth rates for a 10 micron radius droplet are reduced 
to roughly one-half of the growth rates that are calculated when these 
corrections are ignored. 

Finally, the assumption of a one-dimensional isentropic 
flow for moist air is valid only when heat sources and sinks are 
absent and when frictional forces and flow turbulence are ignored, 
The validity of these assumptions depends upon the operating con- 


ditions of the tunnel considered. 
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CHAPTER 3 


THE DROPLET ENSEMBLE MODEL 


3.1 The System 

The closed system (portion of matter) treated by the droplet 
ensemble model consists of two subsystems -- moist air and an ensemble of 
water droplets. Unlike the case of the single droplet model, the equa- 
tions of the ensemble model account for the contributions of the droplets 
to the total energy, momentum, and mass of the system. 

For the ensemble model, as for the single droplet model, we 
consider the matter within a volume contained by the contraction and 
working sections of the tunnel, beginning at the location of the sprayers 
(x = 0.2m). Heat transfer and frictional dissipation at the tunnel walls 
are assumed to be negligible and steady-state assumptions apply. It is 
also assumed that the air flux injected with the water droplets is neg- 
ligible relative to the air flux that would result without the nozzles. 
An estimate of this flux for flow in the NRC tunnel is given in Appendix 
C. The results indicate that its contribution is negligible, 

Although a realistic drop size distribution for the spray is 
described by a range of droplet sizes, the present model is simplified by 
treating an ensemble of N identical spherical droplets per unit volume. 
Each droplet is characterized by a temperature Tor its radius tar and a 
velocity Uo: The droplet velocity is measured in a frame of reference that 
is fixed to the tunnel walls. The droplets are injected parallel to the 
air flow and the droplets are assumed to be uniformly distributed across 
any cross-section of the tunnel. Furthermore, it is convenient to regard 


the droplets as a pseudo-continuum so that an average droplet density 
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(mass of droplets per unit volume of the air-spray mixture) may be 


defined. The droplets are described by the mass density, Pur 


Pre eS 
le) = ate ON 


When droplet evaporation or condensation occurs, the mixing 
ratio of the moist air will not remain constant. Although vapor is 
continuously introduced or depleted by droplet evaporation or condensa- 
tion, it is assumed that the moist air is well-mixed at all times and that 
the moist air can be characterized by a uniform temperature, T, anda 
velocity, U. 

It is emphasized that the "well-mixed vapor" assumption is not 
realistic for tunnels where the time that is characteristic of the droplet 
motion compares to a time characteristic of the mixing, Ta Li hee aks 
assumed that turbulent mixing processes dominate over the molecular 


processes, a characteristic time for the mixing is, 


where Ka is the coefficient of eddy diffusivity and L is a scale or 
characteristic width of the tunnel. An estimate of the time character- 


istic of the motion is given by, 


_ dU 
Te Sy (ee 
For flow in tunnels where an . Oy (i.e. large droplet accelerations), 
adequate mixing cannot occur and the vapor associated with droplet 


evaporation or condensation is not at the same temperature or velocity 


as the surrounding well-mixed air. 
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When the droplets are identical, the droplet growth, thermo- 
dynamic, and drag equations of the single droplet model describe the 
_ Motion and growth of the ensemble. Since the equations describing the 
flow of the moist air and droplets are coupled, the properties of the air 
and droplets are obtained by simultaneously solving the droplet equations 
(2.4.4, 2.5.1, 2.5.7) and the global conservation equations. The global 
equations include two relations that describe the coupling of momentum 
and entropy. In addition, the mass continuity equations for the water 
substance and for the dry air are required. 

The problem of the coupled motion can be considered with respect 
to a control volume fixed in space (Eulerian frame of reference). 
Consider the volume, V, sketched in Figure 3.1. This control volume of 
thickness, 6x, is bounded by the cross-sectional areas A, = A(x) and 
A, = A(x + 6x) and by the lateral surface, Az, which coincides with the 
flow boundary of the system. The steady-state integral equation for the 
mass, energy, and momentum flow through the control volume are determined 
by applying Reynolds Transport Theorem to the matter within the volume 
(Crowe, 1976, p. 387). The theorem states that the rate of accumulation 
of some extensive property (mass, entropy, energy) in a closed system 
is equal to the rate of increase of the property in the control volume 
plus the net influx of the property across the control surface for the 
very instant that the matter occupies the control volume. Mathematically, 
the theorem is written, 


Planes MU Go a (3.1.1) 
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FIGURE 3.1 


: A control volume bounded by the cross- 


sectional areas, A), and Ao, and the 


lateral surface, A3. 
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where a2 is the total time derivative of some extensive quantity, Q, 

and the specific quantity, q = DQ/Dm. The thickness, 6x, is determined 
by the requirement that significant changes in the averaged properties 
should not occur on a scale that is smaller than the length, 6x. The 
scale of these significant changes is considered to be large relative to 
the scale of local changes (i.e. mean free path). Thus, the control 
volume must contain many droplets in order to define an average droplet 
velocity or temperature, for example. 

The mass continuity equation for the water substance in the 
control volume must account for the total water content of the droplet 
ensemble and of the moist air. Since not only steady-state is assumed, 
implying am /at = 0 , but also conservation of water substance is assumed, 


implying Dm, /Dt = 09, 1, tOLlows that, 


fara oon +o Ue)edk .= 0 


In addition, it is assumed that mass transfers do not occur at the lateral 
surface, Az, of the control volume. The steady-state mass continuity 
equation for the water substance that is obtained by applying Reynold's 
Theorem is, 
= + UA wei2 
cata Fe eee} hipaa aa. C 5 (3 ) 
In addition to mass continuity, the number of droplets crossing 
a cross-sectional area must be conserved. Thus, steady-state changes 
in the mass of the liquid water substance occur only by evaporation or 
condensation of a constant number flux of the droplets, The continuity 


equation for the number of droplets is written, 
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N,U_ A, = N,U_A (3.1.3) 
Cc 


3.2 The Global Entropy Equation 

Application of the Second Law of Thermodynamics is simplified 
for a coupled air and water system if specific entropy is used. Specific 
entropy is a well-defined thermodynamic property for each component of 
the system when an equilibrium state exists or when reversible changes 
in its state occur. Iribarne and Godson (1973, p. 75), for example, 
derive an expression for the entropy of a closed two-phase system 
consisting of dry air, liquid water, and its saturated vapor. The 
entropy equation is written in terms of two other thermodynamic or state 
variables, which are taken to be temperature and pressure. 

In Appendix D, an approximation is derived for the entropy of 
a system consisting of dry air, water vapor, and liquid water in the form 
of droplets. The derivation assumes that the temperatures of the moist air 
and droplets are all equal and that the vapor pressure of the moist air 
is the equilibrium vapor pressure defined with respect to a curved droplet 
surface. The resulting entropy equation is an approximation since the 
equilibrium vapor pressure of the droplets is replaced by the saturation 
vapor pressure defined with respect to a plane surface. Calculations for 
a 10 micron droplet at a temperature of 20°C indicate that droplet curva- 
ture increases the equilibrium vapor pressure above the saturation vapor 
pressure by a factor of 1.000086. 

The present system, where the moist air and the droplets move 
with respect to each other, is not characterized by an equilibrium state 
since the temperatures of the moist air and droplets can differ from one 


another and the partial vapor pressure of the moist air can depart 
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significantly from the droplet equilibrium vapor pressure. For such a 
nonequilibrium open system in which irreversible processes are occuring 
(i.e. entropy is not conserved), the laws of classical thermodynamics 
provide a set of inequalities describing the direction of the changes in 
its state. 

If it is assumed that the entropy within a control volume 
remains constant, an entropy equation may be derived that describes the 
state of the coupled air and droplet flow. As in the first model, the 
derivation requires that frictional forces and tunnel turbulence be 
ignored and that the tunnel walls are insulated. The validity of these 
assumptions depends upon the operating conditions of the tunnel considered. 
The constant entropy assumption also means that phase changes between the 
air and the spray are reversible. 

The result of these derivations for the steady-state flow, 
given in Appendix E, is that the total entropy flux crossing a cross- 
sectional area of the tunnel is constant. The constant entropy flux is 


given by, 


a (7) A 
pqUAtc ent ~ R enP 4) + p UA Ee R an ar, 


+ p UAC, ent + p U Ac ent. = constant (32a) 


3.3 The Global Momentum Equation 

When the droplets and moist air have different velocities, 
any phase changes will result in a change of the momentum of the system. 
The integral (area averaged) momentum flux for the steady-state flow of 
the moist air and droplets is determined by applying Newton's Second Law 


of Motion and the Reynolds Transport Theorem to the control volume 
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(see Figure 3.2). 

"The net force acting on the material within the control 
surface is equal to the increase of momentum flux of the streams flowing 
Chrough ther controk surface". (Shapiro, 1953, p. 224). The forces that 
combine se produce the momentum changes are: 

l. the pressure forces at the end cross-sectional areas where such 
pressure differences tend to increase the speed of the moist air; 

2. the distributed pressure forces exerted by the tunnel walls on the 
fluid and producing the required changes in the direction of the 
flow; 

3. the drag force exerted by the droplets on the air. 

Considering the x-component of the pressure forces acting on 
the internal surfaces of the control volume, the net force acting on the 
cross-sectional areas, Ay and Ayr is (PjAy - PA.) « Since the x-compon- 
ent of the distributed pressure forces acting on the bounding surface 
A, is difficult to determine, an average pressure, P = (Cs P.)/2 ne 
used. This approximation is valid for a control volume where the x-com- 
ponent of pressure varies linearly (or more slowly) with the x-distance. 


Referring to Figure 3.2, the distributed pressure force is, 
ieeePydsn @- © f dA =P (AS oA) 


where dA is the projection of an element of the lateral surface on to a 
cross-sectional plane, i is a unit vector in the positive x-direction and 
n is the inward unit vector normal to the lateral surface. The net 


pressure force acting on the matter within the control volume ni 


a : ee - Pp (33300) 
PAA, = PAA, = - (A, A.) A (P, 9) 
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FIGURE 3.2 


: Pressure force acting on a control volume. 


The unit vector in the x-direction is i, and 
n is the unit vector normal to the lateral 


surface, A3. 
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The steady-state drag exerted by the droplets on the air for a 


unit volume of the system is given as, 


N [6rr u(U-u_) (1+are®)] (3713.2) 


F 
drag/volume 


In this study, the momentum transfer of the droplets to the air may be 
neglected relative to the pressure gradient force. Calculations using 
the accelerations obtained from the single droplet model for flow of a 
10 micron droplet in the NRC high speed icing tunnel indicate that the 


ratio of the drag force to the pressure gradient force for 2.4 x 10? 


droplets per cubic meter should be less than igi Thus, the momentum 
changes that are associated with the phase changes are negligible. 
The net one-dimensional momentum flux of the streams flowing 
through the penta) volume may be written as, 
ee + 0 UZA + acs = Pave 7 p UPA - ete 
aan Rp) 
= momentum flux (37373) 
Using this result and neglecting the drag and gravitational terms, the 
momentum equation is, 


2 2 2 Z Does arp 2 
Dae a Oe i 2| a at Wig | A, 
X=X, X=X_ 


= A (P,-P,) (3.3.4) 


CHAPTER 4 


DISCUSSION OF RESULTS FOR NRC HIGH SPEED ICING TUNNEL 


4.1 The Approach 

In this chapter, the results obtained by modeling the flow 
through the contraction and working sections of the NRC high speed 
icing tunnel are reported and discussed. A parametric study is performed 
in order to determine the influence of initial conditions on the flow 
properties at the working section. For the duct geometry of the NRC 
high speed icing tunnel (discussed in Section 1.3), the following duct 
inlet properties are varied: 

1. the moist air thermodynamic and kinetic state variables at the 
plenum (such as pressure, temperature, relative humidity and air 
speed) ; 

2. the liquid water content of the spray at the sprayers; 

3. the water injection velocity and temperature; and 

4. the droplet radius. 

Baseline values of these variables are selected and each of the 
inlet properties is varied one at a time. The following baseline values 


of the initial parameters are used: 


liquid water content of spray 10g/m? 
spray temperature 20°C 
representative radius 25 microns 
spray injection velocity 3 m/s 

air speed at plenum 18,2 m/s 
air temperature at plenum 2022°C 
relative humidity at plenum 50% 
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For an initial or plenum air velocity of 18.2 m/s, a high 
speed air flow is described for the NRC icing tunnel. The downstream 
air velocity at the working section is near 120 m/s. The air temper- 
ature, which is initially 2.22°C, decreases to a value near 1°C at the 
sprayers and has a value of -4.8°C at the working section. The dist- 
ance from the plenum to the sprayers is 0.2 m and the distance from the 
plenum to the working section is 2.1336 m. 

The spray velocity increases rapidly during the first 0.2m 
of the contraction distance from the sprayers and approaches the air 
velocity. The downstream spray velocity at the working section is 
near 106 m/s, indicating that it lags the air velocity by 14 m/s. For 
an illustration of the variation of the temperatures and velocities 
of the spray and air with the contraction distance, the reader is re- 
ferred to the upper curves of Figures 4.11 and 4.12. It is emphasized 
that the plenum velocity in these Figures is 20 m/s rather than 18.2 m/s 
(which is taken as the baseline value). 

The reason for choosing the liquid water content of the spray 
to be 10 g/m? rather than a lower value more realistic of natural icing 
conditions (in the range of 0.5-2 g/m?) is because it is desirable to 
observe the maximum coupling effect of the droplets and the moist air. 
The droplet radius is also chosen larger than a more representative value 
of say, 10 microns. Since larger droplets cool more slowly than smaller 
droplets, the choice of a larger baseline droplet radius means that axial 
variations of other properties are then better observed. The choice of 
a spray temperature of 20°C is determined from the fact that the spray 
water is usually kept at room temperature. The air that is injected 
with the droplets is heated to approximately 30°C to preclude freezing 


of the droplets on or within the sprayers. Although a typical cloud 
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environment is saturated, the baseline relative humidity for this 

icing simulation is taken as 50 percent. If a plenum relative humidity 
is chosen near saturation, the calculated relative humidities rise 

to large supersaturation values near the working section of the tunnel. 
(with relative humidities approaching 150 percent). In a real icing 
tunnel, contamination of the moist air with particles would limit the 
supersaturation. 

The results described in this chapter are obtained under an 
assumption that the droplet disrribution consists of N identical drop- 
lets. A study by Willbanks and Shultz (1973) compared the effect on 
the downstream spray properties of mathematically simulating the water 
spray cloud with a realistic size distribution of droplets to the effect 
of simulating the cloud with a single droplet size. The single droplet 
size was equal to the mass median radius of the distribution. The 
mass median radius is defined as the radius of those droplets for which 
one-half of the entire liquid water content in the spray is contained 
in droplets whose radius is less than this median value. The distrib- 
ution used was a measured distribution determined from a laser holo- 
graphy technique. The size distribution of their study was fairly 
similar to the size distributions found in cumuliform clouds. The re- 
sults from their calculations indicated that the spray temperature was 
not highly dependent on wheth r the spray cloud was characterized by 
having a distribution of droplet sizes or a single uniform droplet size. 
It was concluded that, for typical icing conditions, the character- 
ization of the icing conditions by a liquid water content and a mass 


median radius is sufficiently accurate for testing purposes. 
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The results of the parametric study are discussed in the 
following sections. For a discussion of the details of programming the 
two models, the reader is referred to Appendix H. Values of constants 
and physical relations needed to calculate same of the intermediate var- 
iables used in the dynamical and thermodynamical equations are given in 


Appendix G. 


4.2 The Effect of Variation of the Plenum Relative Humidity 

The effect on the spray temperature of a variation in the 
plenum relative humidity from 0 to 100 percent is shown in Figure 4.1. 

As the Figure indicates, the spray temperature is highly dependent upon 
the tunnel relative humidity. Thus, evaporation and condensation at the 
droplet surface play a major role in determining the spray temperature, 
By controlling the plenum relative humidity, the rate of mass transfers 
and consequently, the spray temperature at the working section may be 
controlled. 

The calculations indicate that the moist air temperature and 
velocity are effectively independent of the tunnel relative humidity. 

The explanation is that the water vapor is only a small part of the mass 
flow of the air. 

Figure 4.2 shows the effect that a variation in the plenum 
relative humidity has on the temperature difference between the spray and 
 Seroptodap It is seen in this Figure and in Figure 4.1 that the spray 
cools at a slower rate for the higher relative humidities since the 
droplet evaporation is less, When the plenum relative humidity is greater 


than 50 percent, Figure 4.2 indicates that the spray temperature lags 
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FIGURE 4.1 : The effect of initial relative humidity on the spray 


temperature as a function of distance from the plenum. 


The arrow indicates the beginning of the working section. 
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FIGURE 4.2 : The effect of tunnel relative humidity on the 


temperature difference (T-T) between the spray 
and air as a function of distance from the plenum, 


The arrow indicates the beginning of the working 


section. 
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the air temperature for the entire contraction distance. For the case 
where the air flow is initially dry at the plenum, the droplet evaporation 
_is greatest and the spray actually becomes cooler than the air flow 
before reaching the working section. 

It is suggested here that the temperature of the droplets 
approaches the wet-bulb temperature as the spray velocity approaches the air 
velocity. The wet bulb temperature, To is defined (Rogers, 1976, p. 21) 
as "the temperature to which air may be cooled by evaporating water into 
it at constant pressure until saturation is reached", From the approximation 


(Ibid, p. 22), 


sere 2 
Be ee Wis ght No 
Save sda: = 
where r is the mixing ratio corresponding to saturation of the air , 


sat 


it is seen that the wet-bulb temperature of the air increases as the mixing 
ratio and temperature of the air increase. When the relative humidity 
of the air is less than 100 percent and r < Coat! it follows that 
RB < T. When the relative humidity is greater than 100 percent (super- 
saturation), the approximation implies that si ars 

Unless the air is saturated, the air and spray (i.e. droplet) 
temperatures will be different. For low relative humidities, it is 
reasonable to expect the drop temperature to become cooler than the air 
temperature since T < T,. When the moist air is saturated at the plenum, 
on the other hand, the tunnel relative humidity rises above 100 percent 
for the tlow along the contraction section. (This effect 1s largely the 
result of the decrease in the air temperature). Since T > Teiaetnis 
case, the droplet temperature can be expected to lag the temperature of 


the air (i.e. the droplet temperature exceeds the air temperature). 
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Hence, the relations shown in Figure 4,2 are explained if the droplet 


temperature tends to approach the wet-bulb temperature of the air. 


4.3 The effect of Variation of the Liquid Water Content 
Figure 4.3 indicates that a variation in the liquid water 
content of the spray has only a small effect on the temperature of the 
droplets. The Figure compares the effect of introducing a large liquid 
water content of 50g/m> to that of having a negligible liquid water 
content. Calculations show that the liquid water content also has a very 
small influence on the temperature and velocity of the air. As shown in 
Figure 4.4, the temperature difference between the spray and the air does 
not vary significantly with the liquid water content. In short, these 
results indicate that there is no significant thermal feedback between the 
spray and the moist air. The reason for this small effect is that only 
a small quantity of water is sprayed into the tunnel relative to the mass 
flow of the air. It should be emphasized that this conclusion applies 
only for the small liquid aha contents typical of icing simulations. 
The solutions given by the single droplet model and by the 
ensemble model are, for practical purposes, identical. The advantage 
of using the ensemble model here if that the CPU time required to run 
the ensemble model on the AMDAHL 470 V/7 at the University of Alberta 
is roughly two-thirds of the time required to run the single droplet 


model. 


4.4 The Effect of Variation of Initial Spray Temperature 
Figure 4.5 indicates that the initial spray temperature has only 
a small effect on the final spray temperature of the working section, It 


is also seen in Figure 4.6 that the final temperature difference between 
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FIGURE 4.3 : The effect of initial spray liquid water content on 
the spray temperature as a function of distance 
from the plenum. The arrow indicates the beginning 


of the working section. 
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The effect of initial spray liquid water content 
on the temperature difference (T-T) between the 
spray and air as a function of distance from the 
plenum. The arrow indicates the beginning of the 


working section. 
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FIGURE 4.5 


The effect of initial spray temperature on the spray 
temperature as a function of distance from the plenum. 


The arrow indicates the beginning of the working section. 
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the spray and the air shows only a small dependence on the spray injection 
temperature. 

Since the calculations from the single droplet model give 
similar results to those from the ensemble model, the earlier conclusion 
that there is no significant thermal feedback between the spray and the 
air is supported. An explanation for the result that two droplets of 
equal size and starting with very different temperatures can end at the 
working section with very similar temperatures is that the droplet 


temperature approaches the wet-bulb temperature of the air. 


4.5 The Effect of Variation of Droplet Radius 

Figures 4.7 and 4.8 show the effect of varying the radius of 
the droplets. As seen in Figure 4.7, a spray of 10 micron droplets 
cools faster than a spray of 25 micron droplets. Figure 4.8 shows that 
the 100 micron radius droplets tend to lag the air temperature by 
substantial amounts for the entire length of the contraction section. 
The smallest droplets, having a 10 micron radius, actually cool below 
the air temperature. Since the calculated relative humidity of the air 
remains below the saturated value for the entire contraction section 
(for a plenum relative humidity of 50 percent) and the wet-bulb temperature 
then remains less than the air temperature, these results suggest that 
the 10 micron droplet approaches the wet-bulb temperature of the air (and 
equilibrium) sooner than the larger droplets. 

The calculations indicate that the larger droplets lag the 
air speed by significantly greater amounts than the smaller droplets. 
The 100 micron droplets, for example, lag the air speed at the working 


section by 38 m/s while the 10 micron droplets lag the air speed by only 


er | aes we ae a 7 
1 Oe ha i i 7 ol 
; y SZ 
; ‘i i 
ee ; Ta 
ri 


aos gnate year at? on cout 18 


ayy fubon setgeiuh aborts a nog “é oad | 

a ; 
reales! yo wektena set! eRe ache ead ee pa ” 
ere? Su. yesge Sit eee! AateTAe ‘scpaeil and 7 
eer we Jad? ten ot a RS eae 
“id a Sine BED BPISISToT ae ee eet yaar mashopionet vale 
Jalgoah. ei Sen? ef srore simnets vst imig’ ie fie 


eC bath ) tensioned aa rs pad ; piety or . - 


ma a olga ta 

dt thytey> nh, POSTS: See “weak i hy ih alae! 

ate gunn sitigin) Of. 3a yang e108 10M ae a adail il 

Suis ignet A cant cele oe a ae 

j ey odt ‘civ ate itt Set - th a m 

an “ee: Oeil eee edt |b ebohashacs ¥ Nic aut? ere f 

ited laos. vitsatnes. ~sthue neo ae shkuae- gman 

“hh AekY Zor GCG ai yee, Geeete ey les! ae eet seracivhcones 

cag ote povikte er Adley ins in is 4 sieldy Bedi eats yet 2 
Sudiotgaed ite air St Ears edi adele Os. 45 oS het aibnetivs 

i ith dean ad Loineo omady nti she “ails ‘Mey deed eater 

) fire) AG GS feoeran bs 196yh30 “pen at dtinseagar a yal is maalh 

cod (hc aaa eats revit ‘aac | 

SY Pl MISE Het SH. cts sist yen a whe 

priigtiow atts aa fig ole coty Fe ol riage ‘en vee 


yin “: eece oie eat, eh 2 altpah vont. ald sae toe 


tei ote Git cab ace nis OMB : : 


2 : J . = ea {| 8) - “a) 


TEMPERATURE (°K ) 
276 280 284 288 292 


272 


268 


264 


FIGURE 4.7 
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The effect of droplet radius on the spray temp- 
erature as a function of distance from the plenum. 
The arrow indicates the beginning of the working 


section. 
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FIGURE 4.8 : The effect of droplet radius on the temperature 


difference (TT) between the spray and air as 
a function of distance fram the plenum. The arrow 
indicates the beginning of the working section. 
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6 m/s. Thus, if it is assumed that the spray reaches the temperature and 
speed of the air at the working section, as is the case in many icing 


simulations, Significant errors could be incurred. 


4.6 The Effect of Variation of the Injection Velocity 

The effect of varying the spray injection velocity is shown in 
Figures 4.9 and 4.10. It is seen that, for these initial conditions, 
varying the injection velocity has only a small effect on the final spray 
temperature. It is noted though, that a smaller initial relative velocity 
between the droplets and the air results in a slower cooling of the spray 
because heat and mass transfers from the droplets are less initially. 
Since the injection velocity has a negligible effect on the air tempera- 
ture, the difference between the spray and air temperatures at the working 
section has only a small dependence on the injection velocity. This 


result is indicated by Figure 4.10. 


4.7 The Effect of Air Speed on the Droplet Properties 

The effect that a variation in the air speed has on the spray 
temperature and velocity is shown in Figures 4,11 and 4.12. The spray 
is injected into the air flow at a speed of 3 m/s in all three cases and 
the slow, medium, and fast tunnel air speeds correspond to plenum speeds 
of 5% 10, and 20 m/s, respectively. Figure 4.11 indicates that the velo- 
city difference between the spray and the air is greatest for the fast 
air flow, with the droplets moving at essentially the air speed in the 
case of the slow air flow. 

For the slow and medium air speeds, it is seen in Figure 4.12 
that the spray temperature rapidly decreases during the first 0.5 mof 


the contraction distance and reaches temperatures that are below the air 
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FIGURE 4.9 : 


The effect of spray injection velocity on the spray tenp- 
erature as a function of distance from the plenum. The 


arrow indicates the beginning of the working section. 
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The effect of spray injection velocity on the 
temperature difference (TT) between the spray 
and air as a function of distance from the plenum. 
The arrow indicates the beginning of the working 
section. 
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The effect of air speed on the one-dimensional 
spray velocity. The air speed is the upper curve 
in each case. The slow, medium, and fast air 
speeds correspond to plenum speeds of 5, 10, and 
20 m/s, respectively. The double arrows indicate 
velocity lags near the working section. 
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FIGURE 4,12 


The effect of air speed on the air and spray 
temperature. The air temperature is the solid 
curve in each case. 


The slow, medium, and fast 
air speeds correspond to plenum speeds of 5,10, 
and 20 m/s, respectively. 


The double arrows 
indicate temperature differences near the working 
section. 
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temperature. In each of these two cases, the droplet temperature varies 
Slowly during the last half of the tunnel length, thus indicating an 
approach to the wet-bulb temperature of the air. When the spray is in- 
jected into the fast air flow, however, its temperature decreases stead- 
ily but is greater than the air temperature for the entire contraction 
length. In the case of the fast air flow, the spray temperature departs 
significantly from equilibrium with the air temperature. 

When evaluating the effects on the behavior of the spray 
temperature of a variation in the tunnel speed, it is difficult to 
isolate the factors that explain the differences. The relations indi- 
cated in Figure 5.12 reflect the interaction of a number of factors that 
determine the state of the air flow and consequently, influence the rate 
of energy and momentum transfers from a spray droplet. For example, two 
main factors that compete to determine the changes in the wet-bulb 
temperature of the air are the increase in the tunnel relative humidity 
with contraction distance, implying an increase in the wet-bulb tempera- 
ture, and the decrease in air temperature with contraction distance, 
implying a decrease in the wet-bulb temperature. The net effect of these 


two factors can have a large influence on the droplet temperature. 


4.8 Summary of the Results 

The following conclusions are drawn by this parametic study. 
The temperatures of the spray and moist air at the working section are 
practically independent of the initial spray temperature, The effect of 
varying the liquid water content of the spray is small, The reason for 
these results is that the spray and vapor fluxes are only a small fraction 


of the total mass flow. Thus, there is no significant thermal feedback 
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between the spray and the air. The droplet temperature is largely 
controlled by properties of the surrounding air and tends to approach the 
wet-bulb temperature of the air. It is found that a substantial control 
over the thermal state of the spray at the working section is possible 
through control of the plenum relative humidity, A large initial rela- 
tive velocity between the injected spray and the air helps to promote 

a more rapid approach of the spray to thermal equilibrium with the air. 
In addition, droplets having a radius of say, 10 microns, approach the 
wet-bulb temperature of the air at a much faster rate than larger drop- 
lets having a radius of 100 microns. Thus, a shorter duct length can be 
used for the smaller droplets to achieve an effective equilibrium with 
the air. Finally, for plenum air speeds of 20 m/s, the temperature and 
the velocity of the spray significantly lag that of the air. For plenum 
air speeds of 5 and 10 m/s, the spray approaches the velocity and wet- 
bulb temperature of the air, while it remains substantially out of 


equilibrium with the air in the case where the plenum air speed is 20 m/s. 
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CHAPTER 5 


CONCLUSIONS AND RECOMMENDATTONS 


5.1 Final Conclusions 

The present study has developed a model of general applicability 
for determining the actual spray and air dynamics and thermodynamics in 
any testing facility. MKS units were used for all variables and constants. 
The properties of the spray and air, averaged over a cross-sectional 
area perpendicular to the tunnel axis, were determined using a one- 
dimensional numerical model of the flow. The main properties calculated 
were the temperatures and velocities of the spray and air, the droplet 
size, and the air pressure. The single droplet model and the ensemble 
model were developed for this purpose and the governing equations were 
programmed in FORTRAN IV language for solution on a digital computer. 

The single droplet model, introduced in Chapter 3, applies to situations 
where the spray has only a negligible influence on the air flow. The 
ensemble model, discussed in Chapter 3, is a modified version of the 
simpler single droplet model, applicable to cases where the contributions 
of the spray to the moist air flow have to be considered. 

A limited parametric study was performed to determine the 
effects of the initial moist air and water spray conditions on the state 
of the flow at the working section of the NrRc_ high speed icing tunnel. 
A number of significant results were obtained from this study. It was 
found that the final solutions given by the two models were essentially 
identical for the liquid water contents typical of an icing situation. 

It was concluded that there was no significant thermal feedback between 


the spray and the moist air. The spray temperature at the working section 
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of the tunnel was found to be sensitive to the relative humidity of the 


hb ae 


It was considered that, for sufficient duct length, the droplets 


tended to approach the wet-bulb temperature of the air, with 10 micron 


droplets approaching this temperature more quickly than the 100 micron 


droplets. For plenum air speeds of 5 and 10 m/s, it was found that the 


velocity of the 25 micron droplets tended to approach the velocity of 


the air and the droplet temperature tended to approach the wet-bulb temp- 


erature before reaching the working section. When the plenum speed was 


20 m/s, comparison with the results for the slower tunnel speeds indic- 


ated that the droplet properties remained significantly out of equil- 


ibrium with the air properties for the entire contraction section. 


Sc 


Final Reconmendations 


The following recommendations are given, some of which have 


been mentioned in previous discussions: 


te 


Droplet freeze-out processes should be incorporated into the droplet 
growth and thermodynamic equations. In this respect, the inclusion 
of particulate matter with the moist air flow would contribute to 
simulation of freeze-out processes and could ensure that more 
realistic tunnel supersaturations result when the plenum relative 
humidity is near saturation. 

Calculations performed for flow of a single droplet in the NRC high 
speed icing tunnel indicate that the history term should be included 
in the droplet momentum equation when the injection velocity of the 
droplet relative to the air is small. This result was obtained 
using the droplet accelerations that were calculated by omitting 


the history term from the equations of motion. 
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Consideration should be given to relaxing the assumption that the 
vapor introduced by evaporation of the spray is immediately mixed 
with the ambient moist air. It is suggested that an entropy flux 
equation might be derived that accounted for two vapor tempera- 
tures —- one temperature associated with "recent" changes in the 
vapor density introduced by evaporation and condensation of the 
spray and one associated with the well-mixed portion of the moist 
air, 

The present ensemble model should be expanded to incorporate a non- 


uniform size distribution. 
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APPENDIX A 


LIST OF SYMBOLS 


flow cross-sectional area 
specific heat capacities at constant pressure for 


dry air and water vapor, respectively 


specific heat capacity at constant volume for dry 


air and water vapor, respectively 


specific heat capacity for water 

emissivity of water 

partial vapor pressure of moist air 

saturation vapor pressure over a droplet 

saturation vapor pressure with respect to plane water 


surface 


entropy flux crossing a cross-sectional area 
mean ventilation coefficient 
acceleration of gravity 

thermal conductivity of moist air 
specific latent heat of evaporation 
mass 

mass of a droplet 

mass of system 

Mach number 

number of droplets per unit volume 
mean Nusselt number 


Prandtl number 
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pressure of moist air 

mass flux of dry air and water substance, respectively 
specific humidity 

Specific gas constants for dry air and water vapor, 
respectively 

Reynolds number 

relative humidity 

mixing ratio 

droplet radius 

entropy of spray and moist air system 

Schmidt number 

mean Sherwood number 

specific entropy of dry air, vapor, and spray, respectively 
temperature of moist (and dry) air 

virtual temperature 

time 

one-dimensional velocity (moist air) 

one-dimensional spray velocity 

axial distance from plenum 

mass flux of air injected by a sprayer nozzle 

dynamic viscosity of air 

kinematic viscosity of air 

mass density 

density of water 

ratio of specific heat capacities for moist air = Coy 
Stefan-Boltzmann constant 


surface tension of water in air 
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APPENDIX B 


MACH NUMBER FORMULATION FOR FRICTIONLESS, ADIABATIC FLOW 


The equations that determine the adiabatic frictionless one- 
dimensional flow of dry air through the contraction of the wind tunnel 
can be rewritten in terms of the local Mach number. This non-dimensional 
variable is defined as the ratio of the local flow speed to the speed of 
sound. For adiabatic flow of dry air treated as a perfect gas, the Mach 
number is given from Equation 2.2.8 as, M= B/ RAT) 

The Bernoulli equation for compressible adiabatic flow can be 
integrated to yield Equation 2.2.5, 
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Using the Perfect Gas Law, and the definition of Mach number in the form, 
iF = MoyRT, Equation 2.2.5 can then be written in terms of the tempera- 
ture and Mach number as, 


es T 2 R 
Mo ahalo Vala ee Cn 
na SE OI ee Fe peter ewes nen Oy ae : 
2 Y-1 2 y-1 1 
It is assumed that y is a constant in the above expression. Equation B.1 


may be rearranged to give the relation, 


2 2 & (B.2) 
2 cmiaeibantainaoa + + esa ° 


Substitution of the Perfect Gas Law and the definition of the 
Mach number into the continuity equation, (2.2.4), results in an expression 
relating the temperature to the air density for two points along the flow 
direction. The expression is, 
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Equation 2.2.3 (Poisson's equation), describing an adiabatic 
expansion or compression of a perfect gas, can be written in terms of 
density and temperature by applying the Perfect Gas Law. The result is, 


tT - | y | (yal) (B. 4) 
eva!) ay) C 


Combining Equations B.2, B.3,and B,4, a relation is obtained 
that links the Mach number to the ratio of the cross-sectional areas at 


positions xy and Xo Using y = 7/5 = 1.4, the relation is written as, 
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APPENDIX C 


CALCULATION OF THE SPRAYER AIR FLUX FOR THE NRC TUNNEL 


The object of the following discussion is to show that the 
mass flux of the air that is injected through the nozzles is small 
compared to the mass flux without the sprayers, 

The air that is injected with the spray is contained within 
the spray nozzle between an outer tube having a diameter of 3.1 mm and 
an inner tube having a diameter of 1.3 mm. The water that is injected 
is contained within this inner tube and a smaller concentric (hollow) 
tube having a radius of 0.51 mm. 

The assumptions that are used to obtain an estimate of this 
flux include the following: 

1. The air is flowing at the sonic speed through the nozzle (that is, 
the nozzle is a critical nozzle). 

2. The air is heated to 30°C within the nozzle. 

3. Upon injection, the air expands adiabatically to achieve a tempera- 
of -40°C and a pressure equal to the (controlled) ambient tunnel 


pressure (Lozowski, personal communication). 


A relation may be obtained connecting the mass flux to the 
temperature and pressure of the injected air. The expression for the 


mass flux is written as, 


se (Oxi) 
eal (ae) a 


where c is the sonic speed. The speed of sound for a perfect gas is 
determined from the air temperature and is given by, e7 = yRT, Substi- 


tuting this expression into the mass flux evaluation, the result is, 
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where the subscript n describes variables within the nozzle, The air 
pressure within the nozzle is calculated using the adiabatic expansion 


relation (Poisson's equation), 
Bre Mer P iets 
Sate et ona Cas 


where the subscript i describes variables after injection. The injected 
air flux is determined using Equations C,2 and C.3. 

For the initial conditions given in Section 4.1, the calculated 
tunnel pressure at the sprayers (P.) is 99,710 Pascals. The pressure 
within the sprayer nozzle, given by Equation C.3, is approximately 250,130 


one the injected mass flux is 


Pascals. Using T, = 30°C and A, = 6x10 
obtained from Equation C.2 and is 6x10 °kg/s. The results from calcula- 
tions for flow in the NRC tunnel,on the other hand, give a tunnel air 
flux at ihe sprayers of 13.4 kg/s. The ratio of the injected flux to the 
tunnel air flux, 4.5x107*, indicates that the injected flux may be 


neglected relative to the air flux that would result without the nozzles. 


1 ;Vat 7 
i [ws =m be 
/ io, ae A * 
coy 7 - 
; - 4 7 ie 
i 


c&.. 0) 


sti att aiserok oes wari lw apices 
MO LAKEYYo atyedetie an) rng eet 


fats eat nei caatir “tavae ieee = ee ore . 
go Mars on ro Saw 

7 iz. 7 + 

) Farce “as ‘ace RANT Tey: at 

= wis . 

Walotep eit . ib wobdaee Ne cee Sh ver | 
O20 “edenniwians.ef (fv pines vc Pekaags stones 
7 oy. 
fe ide riba ate Sang , Sg com 
“Si aRLN® (hire boar anr er Tae latin 
16 dstined cody . bed bee sed ee, 
vy : at 

ais cue world) ie si bie +o mh) ingss Shh Re a 
. v : a % . | | 


aon ie Wns 4 PoTiavel “ap leah 4 
oy ? 
| ey 


ie 


APPENDIX D 


EQUILIBRIUM ENTROPY FLUX 


According to the Second Law of Thermodynamics, the specific 
entropy of a system in equilibrium is a well-defined function of its 
state. Given a system in equilibrium consisting of Ms kilograms of dry 
air; mM, kilograms of condensed water in the form of droplets and mh 
kilograms of the saturated vapor, the specific entropy of each of the 
above components is a function of two independent thermodynamic (or state) 
variables. Accordingly, expressions can then be written for the contri- 
bution of each component to the total entropy. Letting Sr denote the 
total entropy of this system, letting Sar So and So denote the entropy 
of the dry air, vapor and droplets, respectively, and letting Sar Sur and 
S. denote the corresponding specific entropy quantities, the total 


entropy of the system is given by, 


So =mese ms. ms (D.1) 


For an equilibrium (reversible) change in this state, the change in the 


total entropy is written, 


= + 
DS DS | + Ee DS. 


where D refers to the total differential. For each component, 


DS 4 = msDS 4 

ps0 =e nieDs ee Sem (D.2) 
V V 

DS) = m-D0s. + s. Dm 
(@. Cc 


78 


6 


3 Ace 7 i ty se 5 a : 
me i td a Na 
1 mn ieee Mt fas 
et a se 4 or ie if 
bert 
wots Fy aera i wea ; iat 
ae ha Wy vi 
ar) i ; LD : 7 wy - : : ; a 
Pom i) a ; 
vate Seine te sin) babaee oe w ee vet 1 ae 
7 7 a “i “a 
4 ce +e c oes Ba) Ci a ek es ‘a Fe Ste rie OP p of Dea 
- f Ph F ie rh a “— 
Ma i epi } T Ligh ie Fates J rir}. « 
7 . ‘ oe 
‘ a ma? ray bey Pr. seater oO nT ete | 
¥ vara t r et eS SH hee tat Ke » i - ey |: a Pd 
i wyney ale Le fit , aoe ies ge tnaine ants bs § 
4 a 7 ! : . I 
5 t wa ay fi, ad Kare, 6 Si i baht HK, 
‘c rte BD eet] eaieealll rete dame tod al att 
/ ' Mi p ee i: 7 “a. 
tt NIL Verne: LT ee 14 dk wenn Mm F tan a) a 
: re On 
7 2 ge OES yy RP ay ; vat Peak fe ee alate aa) ee reins ae 
nd v ‘ | Le ( | ey 
ee 8 isl mis Vievieiegie: s7Elen ‘i? Sir Tite veh a 9 
: : t ‘ey. ohn @ Pr a L 
: re 
ms Lv bal BI + } aie VRROSS 4 ars <it te ign: “er , at ; #200 b ie 
M al ao aa 
whieh ti aufaitt age eet 
‘ [re 
a p j + i" 
: (ica boa b M 
) een er, 
ee ll aT 
ae 
; mvs oe - _ adi 7 
Ay ae Vosgesayetss ao Sipe. bites are Iv Piltathe me Ot 
i 7 if P ve 4 th ‘im i 
¢ “i aie - 1 Youn ¥ P nace 
em ial we! yatnetaa “lags 
2 “—_ ‘ a ae ea 
, aoe a + ee rr 
“4 » i , Pd a i _ i: 


19 
where Dm = -Dm.,. If pressure and temperature are the state variables 
used, it is said that an equilibrium expansion or contraction of the 
system occurs if the temperatures of all three camponents are equal to 
each other and if the vapor pressure of the air is maintained at the equi- 
librium value for the droplets. In order to simplify the derivations, 
the equilibrium vapor pressure for the droplets is approximated by the 
saturation vapor pressure with respect to a plane water surface. Calcula- 
tions for a 10 micron droplet at a temperature of 20°C, for example, 
indicate that the factor by which curvature increases the equilibrium vapor 
pressure above the saturation vapor pressure is 1.000086. 

According to Iribarne and Godson (1973, p. 75), total differen- 
tial expressions are given for the change in the specific entropy of 


each component if changes in the variables occur under equilibrium 


processes. The relations are, 


a5 aS | 
ga ut dP = D(gn T) - R,D(en P 
Ds | 7p 7 dT + ai ie P oe (one) 4 (Qn a 
aoe | ae AP 
J ee 
Spee aban dP = D(gn T) - R_D(gn P 
Ds P|. dT + ar bil Ge oy (gn T) (Qn 4 
ar fia» 
os 
jor | dT = c D(2n T) (D.3) 
. 9 Pim 
a 
as ¥ 
where (Spo pm * 0. 


The subscript Tp is used to describe partial differentiation with the mass 
composition of the system held constant, 

It can be shown that the specific heat quane es iC and C azt 
remain constant if the gas behaves according to the Perfect Gas Law 


(Shapiro, 1953, p. 42). If it is also assumed that the specific heat of 
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water remains constant for the small change, D(%nT), Equations D.3 are 


easily integrated to obtain the specific entropy of each component. The 


result is, 
Pp 
T d 
Sy = Sy +c Pilg R erp 
0 0 d 
0 
(T) 
§ .= 5s +c pene Ren =e - 
eg Our rae 0 sat’ 0 
So = 5s + cone 
= 0 (D.4) 


The subscript 9 denotes a quantity measured (or calculated) at some 
reference equilibrium state. 

Using Equations D.3, it can be shown that changes in the 
total entropy of the moist air and spray system that is in equilibrium 
faced only on changes in the state variables, T and P,, and in changes 
of composition. In other words, changes in the total entropy are 
independent of the reference state chosen. 

It is useful to introduce two additional relations, the Clausius- 


Clapeyron Equation (Iribarne and Godson, 1973, p. 62), 


sat RS Vv (D.5) 


re ce = c )dT (D.6) 


The quantity, doe in the above relations refers to the specific latent 


heat of vaporization, Using Equation D,.6, the Clausius-Clapeyron 
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equation is integrated by parts to yield an expression for the satura- 


tion vapor pressure. The result of the integration is 


(T) T; Q i 
oo: ae 7 a | ee ee 
ee. Ty - DT ut + te DY: 
0 0 
= sas - sal + (Gig e= co .en Z (D7) 
0 coat 0 
From Kirchhoff's Law in integrated form, 
Q = = 2 
7 oy c) (T Ty) +k 
0 
and Equation D.7, it can be shown that, 
SAN 
Re Oe (D.8) 
s,(T) ap l,) 7 
If the following mass conservations relations, 
m =—COnStant = m + om 
Wo Vv e 
a = constant = Me 
0 
are substituted into Equation D.1, it follows that, 
ae = Di. 
S, m4Sq + m (s s_) +o Se (D.9) 


0 


Using the appropriate specific entropy relations (Equations D.4) and 


Substituting Equation D.8 into (p29); «the result 1s, 


S = = = 
T m4 Ebs er Ty + Ryan i ee Th 
“ Qn T meee) 
m ic en + m, oat t= Ryan Py + roy eae (D..10) 


xe + Ql wy a 


v 


— doit oo ak Ds ofl ¢ 


tf 
. Spospiler, eae sem i oh 
a if | 
a # delves 64 a 
Ny wv iw ail Cele: 
eck 
tes Lene 
(e.0) 


* le Ag An 


; athh ee 


7 7 ie 
‘ i on 7 
a . ac st a, a 
) Lo , vie tia, i a i ul ‘ ; 
7 iad if - f iP | ' ; 
7 tad oo eee ; > a “ ia 
Sey & boa 


82 
The constant terms on the right hand side may be combined to give the 
reference entropy, Sos The following relation results: 
m2 (T) 


Vv 
i acta myR yen Py (Det) 


Equation D.11 states that the total entropy depends only on the state of 


the system and on the composition of the droplet and moist air system, 
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APPENDIX E 


DERIVATION OF AN ENTROPY FLUX EQUATION 


Let us consider the case where the water droplets and the 
moist air are moving in a tunnel of varying cross-sectional area, The 
droplets and the moist air are each characterized by different tempera- 
tures and axial velocities so that a nonequilibrium flow is described. 

In addition, we assume that the partial vapor pressure at any cross- 
section may depart from its equilibrium value. 

The rate of accumulation of total entropy within a control 
volume may be determined by applying the Reynolds Transport Theorem to 
Ene, eneropy tlow. ‘The result ys 
—— = < eS korcg dai etane Reels Ucdh (E21) 

t v lat S Ie ALi aU 
The summation accounts for contributions from the dry air, water vapor, 
and droplets. Since steady-state conditions are assumed, the entropy 
does not accumulate because of internal sources within the control 
volume. Since it is also assumed that the control volume is adiabatically 
insulated at the lateral surface, entropy can not be exchanged at this 
surface (also, the flow is one-dimensional). For steady-state flow in 
the tunnel, the total entropy in the control volume remains constant only 
if no net flux of entropy crosses through the cross-sectional flow areas, 
If it is assumed that the total entropy within the control volume is 


constant, an entropy equation can be derived to describe the coupled state 


of the system. 
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Referring to the control volume in Figure 3.1, the total 


entropy flux that crosses the area, A = A(x), is constant and is written, 


Fe = constant = F. (x) = F(x+dx) (E22) 


where the variable Fo refers to an entropy flux and Fo is its conserved 
0 
value. Using the appropriate one-dimensional velocities, the entropy 


flux crossing the area, A, is given by, 


(E53) 
Fg (x) i (s ye g’ P Ai . ve x 
The specific entropy variables, Sar Sor and S, are defined in Appendix D. 


If the mass continuity equations, 


0 ,UA = Q, = constant 
d d (E.4) 


(o U a: pU)A = Cue constant 


are substituted into (E.3), the entropy flux at any position x is given by, 


F. (x) = em s(T) = Sedde) see oN (E.5) 


The specific entropy variables are, of course, evaluated at the appropriate 
temperatures. | 

Let us define an equilibrium reference state in which the 
temperatures of the droplets and the moist air are equal. We consider that 
the reference temperature, Tor is equal to the saturation temperature 
where, for practical purposes, saturation is defined with respect to a 


plane water surface. The specific entropy quantities at the reference 


V 
0 0 


entropy quantities given in Appendix D, the entropy flux crossing a cross- 


state are denoted as hi S= vand So Using the appropriate specific 


sectional area is written, 
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Se ieaeao! 0,36 (T)] ie legat Tae anes, 
ue (E.6) 
a - E.0 
+ Qc gn i s 0 uals (T) 5. (T.)] . 

The entropy difference, S_(T) - S(T) , for nonequilibrium 
temperatures, T and Tor and for the nonequilibrium vapor pressure, e, is 
given by, 

T) - = - re = 
EP =A y129 (Tp)0>) 2) pl | tevergen E ag-h hte 


Cc 
0 0 


Cc 
= cen pe (E.7) 
The nonequilibrium vapor pressure term and the non equilibrium tempera- 


ture term in the above equations may be separated as follows, 


e e sat 
R&n— = R Ln + R 2&n 
V ey V esa) V eo 
a i 
SS —_>—sv—s+ a 
CG kn Ts, ee T Ts 


Using these results in (E.7), the result is 


- 
re = - + Cee ae 
s(T) S/40) s (TQ) s_(T)) (coy =) as 
0 
e (T) ik 
sat Cc € 
2 eka - ¢c tn =—- R &n 
Ron eo W aa Vv an) 


(E.8) 


Since the reference state is also an equilibrium state, the first term 
is evaluated as, 
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Substituting Equations E.7 and E,9 into Equation E,8, the result is, 


a (7) ir 
o AT) - s_(T,) jae - R In (RH) - cin = (E.10) 


where RH = e/e (T) 


sat 
If Equation E.10 is used in the entropy flux equation, (E.6), 


the resulting equation is written, 


F = )+ - 
5 (x) UyBalultel Spal a) Qye4gen its QR en ia 
- Q cen Ty + Q, fe, gen T- R yen ay + Qc en ifs 
& (T) the 
+ p. UA + = Ran (RH) ~ p,UAc £n 7 a 


The first five terms on the right hand side of ( E.1]) are constants and 
they may be combined to give the (conserved) reference state entropy flux, 


Pe weercuacion 1). then becomes, 
0 


( 
f(x) carbo eed lc In T - RJIn 2 + p ual hy 


\ 
pd d + - AB SEUnE 


+ pYAc In T + pU Ac In lige CERT?) 
For the case of no droplets (o, PAO), sls eneropy flux relation reduces 
to Equation 2.2.3, which describes the isentropic expansion or contraction 
of the moist air. 
In summary, the entropy equation (E.12) states that, for steady- 
State assumptions, the entropy flux crossing a cross-sectional area of an 
adiabatically insulated spray and moist air system is constant. The state 


of the system is then described for all cross-sections and at all times. 
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APPENDIX F 


FINITE-DIFFERENCES AND NUMERICAL PARAMETERS 


Since the systems of equations describing single droplet and 
droplet ensemble motion in an icing tunnel are closed, the differential 
equations may be solved numerically. In this study, the equations are 
programmed in FORTRAN IV language for solution on a digital computer. 
The derivatives are approximated using an explicit Eulerian finite- 
difference scheme. Boundary conditions are specified at the plenum 
(beginning of contraction section) and at the position of the sprayers. 
For a listing of the programs, the reader is referred to Appendix H. 

The explicit Eulerian finite-difference scheme is first order 
accurate and is the simplest finite-differencing scheme available. 

Since the differential equations used in this study are complex, it is 
desirable to use a simple finite-difference scheme, The disadvantage of 
using the Eulerian scheme is that a small grid spacing must be used in 
order to reduce the effects of numerical truncation error and to reduce 
computational instability. The use of a grid spacing smaller than some 
optimum size, however, is uneconomical. The optimum grid spacing, defined 
as the maximum grid spacing needed to obtain solutions of given accuracy, 
is not constant over the tunnel length since it is dependent upon the 
behaviour of the solutions. 

In the present study, the grid spacing is a function of the 
tunnel geometry. Consequently, any scheme for determining grid spacings 
is tunnel specific. In most cases, however, the grid spacing near the 
sprayers has to be very small relative to the spacing at other positions 


further along the axis since large accelerations of the droplets occur 


87 


yn tev ogee vol etl cree siaes BOL Seppe Edt 
epryeSh FD cas . 1 naal > omar. he iotentast gtk rt. 
eve erobmuupe oft ates oer ie sec eg li ek te 
pectoroner rink et iS Roe ACI tae 7% Ayeutaies, vi Sane 
eek geile Pear toe ines nean atest srk: 
nirete 2 te felipe Sher er gsebeneat spear 
dionvuie aft a0 unt tiaty atthe af hale ington’ ndlsoictce te + nad rt 
# «thee ao bower: 2. wae feet panenietelaa 
Yoittes 22 a evevibe sonenettibssd a7 sadvoline steels air” 
ene yt RITES el int (ie bayite owt at wai 
+i) .eiae och yhude eek ley ~aertanar a 
in cipaditey bec th od ated © brides eer aplegeytat: ‘aot ob pt 

rir’ Sed Fe. Tan trap 1G hae T Taunt ieee iy onal anisole it ped 


ee ai . vam 
enue geht 7) Tans ge ogee, i ny |S ons ar ruc | 
Feautes srtieve GG. eon svi SEM, wana) gt Sashes pea: 
parses Toate: hee ta at 
ats dies satan biases BT ey snl thigney” ore vont doug dation: 2 
Sant! sea tan 2 

Mi? 26 Attowst 5 2. enisade biwe get a Fi im be L 
are cance Art pila erecta act cantante. 17 pi se 
689 1 pritopcye Bche A, eee asi seem re eae i 
eral shetywe: ts, witli cium atts od ovens Lema sii et ot aert 
BUD jateaboeats sft 6 oF sr afesos spiel Sorts ry a9 


e 
7+ 
+ 
~ 


Pe | 
. 
m 

' 

My 


rey Oh ae Sta iS gi ae mee 50 Ae 


a\ Neos Tele 2 Biel fy (om ft AS oh : 


A 


a - “” ‘ ® vw 


88 
near the sprayers. 


The vertical cross-section of the contraction section of the 
NRC tunnel consists of a rapidly converging bell-mouth (circular cross- 
section) followed by a more gentle linear contraction. The cross-sectional 


area of the bell-mouth contraction for x < 0.1524 is given by, 


A = 4{0.3810 - (0.02323 - (x - 0.1524)2)?]? 


The cross-sectional area of the linear contraction is described by, 


A = 40.2286 - (0.03846(x - 0.1524)] * 

For both the single droplet and the ensemble models, optimum 
grid spacings are determined from the requirement that the uncertainty 
in the solutions resulting from a choice of the grid spacings is less than 
the uncertainty in the solutions resulting from errors in the measure- 
ment of the tunnel dimensions. The approach adopted is tO increase the 
grid spacings. for successive runs using the same initial conditions until 
the difference between the calculated properties at the working section 
and some "standard" properties exceeds an acceptable error (i.e. the un- 
certainty associated with measurement errors). The properiees considered 
here are the temperature and velocity of the spray and the moist air. 

The "standard" values of these properties are determined from a run that 
uses very small grid spacings and is restricted by a very small error 
limit (to be discussed later). The initial conditions used in calcul- 
ating these properties are discussed in Section 4.1. 

The grid spacing scheme used for the NRC high speed icing tunnel 
employs logarithmic increments for sections of the tunnel near the spray- 
ers and equispaced increments for sections of the tunnel where the per- 
imeter of a cross-section varies linearly with the x-distance. The details 


of this incrementation scheme are given in Table F.1. 
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TABLE F.1 


GRID SPACING PARAMETERS USED TO SOLVE FOR FLOW 


IN THE NRC HIGH SPEED ICING TUNNEL 


X-Range Increment Number of Equation for Incremented 
Type Increments Position 
(m) (m) 
(tnx, + 1077) 
Ooi 202 Log 1000 X, =e 1 
-4 
yh Sage Log 1500 et tl 
Binie-62. 1336 Linear 4000 x, =x, + 46x10" 


The grid spacings shown in Table F.1 are obtained using an 
assumption that the dimensions of the NRC tunnel can be measured to a 
tenth of a millimetre. The results of two numerical experiments that 
vary the tunnel dimensions by amounts of the uncertainty in the measure- 
ments are summarized in Table F.2. (The relative changes given in this 
Table are percentage diviations from the "standard" properties that are 
obtained for the "actual" tunnel dimensions). In the first experiment, 
the radius of the bell-mouth contraction is decreased by .000l mtoa 
value of .1523 m. The second experiment reduces the slope of the linear 
contraction by an amount of .0001 (units of vertical distance per unit 
of x-distance) so that the altered cross-sectional area varies more slowly 


than the actual cross-sectional area. 
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TABLE F,2 


RELATIVE UNCERTAINTIES IN THE SOLUTIONS TO THE ENSEMBLE MODEL 
RESULTING FROM HYPOTHETICAL MEASUREMENT ERRORS. THESE RESULTS 


ARE OBTAINED FOR FLOW IN THE NRC HIGH SPEED ICING TUNNEL. 


Percentage variation in solutions at the working 


Hypothetical section from the standard 
Measurement 
BELO U U iy a 
c eS 
(m/s) (m/s) (°K) (°K) 


Slope of linear 
contraction 
from .03846 to 


03845 erie 4 


3 2 


4.0x10- Br Oscon2 2.0x10- 
Radius of bell- 

mouth contraction 

decreased from a 2 

2 1524M c0, Gb523M Lye BeO%k0 23 Deeks 


The conclusion from Table F.2 is that a smaller uncertainty 
(in the final solutions) results from an error in measuring the slope 
than from an error in measuring the radius. The experiments indicate that 
the temperatures are more sensitive to changes in these tunnel dimensions 
than are the velocities. In order to determine the optimum grid spacings 
with the best precision, the acceptable error for this particular study 
is chosen as the smallest uncertainty appearing in Table F.2, Hence, the 
grid spacings of the NRC tunnel that are given in Table 4.1 are obtained 
from the requirement that the calculated spray velocity at the working 


section differs from the "standard" value by not more than ,0004 percent, 
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oY 
This requirement is more stringent than necessary because of the uncer- 
tainty associated with using averaging techniques. 

In the single droplet model, the equations for the moist air 
flow can be integrated to give a velocity equation, (2.2.10). For the 
ensemble models, on the other hand, the equations describing the coupled 
air flow are written in finite-difference form and are solved simultan- 
eously by iterative techniques for the moist air variables P, T, and U, 

The initial estimate of the velocity needed to solve this 
system of equations assumes incompressible flow. Letting the subscript 
2 denote unsolved variables at the position X5 and letting the subscript 


1 denote solved variables at the position x,, the following steps are 


VY 
used to solve for the unknown variables in the ensemble model: 


1. The droplet properties at position x, are obtained from the droplet 


2 
flow equations (2.4.4, 2.5.1, 2.5.7). These finite-difference equa- 
tions are written in terms of the known moist air and droplet proper- 
ties at the position, Xj. 


2. These calculated droplet properties are used to determine the mixing 


ratio of the moist air at position X5 from the mass conservation 


relations, 
Q, = p UA = constant = %, 
= + Ao = t = Q 
Q| p, UA pave constan Wo 


Since p UA =rQa- the mixing ratio may be obtained using the expres- 
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An initial estimate of the moist air velocity, U5, is calculated from 
the assumption that the air flow is incompressible, 

Using the estimated velocity Une the moist air pressure PS at position 
X5 is obtained by solving the global momentum equation (3.3.4). The 
result is the estimated pressure, P,. 
The estimated variables, U5 and Po, are then substituted into the 
entropy flux equation (3.2.1) in order to solve for the variable, T,. 
An improved estimate of the velocity Us, is obtained by substituting 
PS and TS into a relation obtained by combining the continuity equation 


for the dry air and the equation of state. The relation is, 


QR T*(l+r) 


22 
If the absolute value of the relative difference between the improved 
estimate of the velocity and the preceding estimate is greater than 
an acceptable relative error, EPS2, steps 4 to 6 are repeated until 
the relative difference is less than EPS2, 


For both of the computer models, it is necessary to determine 


a value for the acceptable error limit that terminates the iterations 
needed to solve for the air velocity. In the case of the single droplet 
model, the iterative equation (2.2.10) is considered to be solved when 
successive improved velocity estimates differ by amounts less than EPS2. 
When the equations for the air and spray motions are coupled (i,e. the 
ensemble model), it is seen in step 7 above that the parameter EPS2 plays 
a Similar role. 
The procedure used in determining EPS2 is similar to that used 


in determining the grid spacings. The value of EPS2 is increased for 
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Successive runs until the last increase produces calculated properties 
at the working section that differ fram the "standard" properties by 
amounts exceeding the uncertainty of the solutions, Values of EPS2, as 
well as the range of the absolute errors incurred in iterating for the 
velocity, are summarized in Table F,3. In this study, the value of EPS2 


used is lose 


TABLE F.3 


ACCEPTABLE UNCERTAINTY LIMITS AND THE ACTUAL RESIDUALS 
THAT RESULT FROM ITERATION FOR THE AIR VELOCITY IN THE 


NRC HIGH SPEED ICING TUNNEL 


Acceptable relative Equations Number of Range of Magnitudes 

uncertainty in air solved iterations of the relative 
velocity (EPS2) (n) to solve errors (E) resulting 

for velocity from iteration 

E appl Use| 
U 
n 
ORS (222916) 1°36 hero +) Seto 
Hae Global ae is Dies ae 
Equations 


A test of the accuracy of the solutions to the finite-difference 
equations for the moist air flow may be performed by comparing the results 
obtained from the single droplet model to those obtained from the ensemble 
model if it is assumed that the initial liquid water content is zero. 

Since the moist air equations may be integrated analytically for the single 


droplet model, this comparison provides a reasonably good check on the 
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94 
behaviour of the solutions to the ensemble model. 
For flow in the NRC high speed icing tunnel, calculations } 
from the two models (using identical initial conditions) indicate that the 
two sets of solutions vary by negligible amounts (less than Tome percent) . 


The two sets of solutions are compared in Table F.4, 


TABLE F.4 


A COMPARISON OF RESULTS FROM THE SINGLE DROPLET 
MODEL TO THOSE FROM THE ENSEMBLE MODEL (LWC = 0) 
FOR SIMILAR INITIAL CONDITIONS. THESE RESULTS 


DESCRIBE FLOW IN THE NRC HIGH SPEED ICING TUNNEL 


Model Solutions at working section 
U U T ats 
( e) 
(m/s) (m/s) (°K) (°R) 
Ensemble 119.9976 LOSRDVaS 206.2351 269.9343 
(Numerical Solution) 
Single Droplet 119.9976 TO5.57 72 268.2349 269.9342 
(Analytical Solution) 
Absolute Difference 0 0.0001 0.0002 0.0001 
5 5 5 


Percentage Difference 0 9.5x10— Peaxtoe Sypelels 
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APPENDIX G 


PHYSICAL CONSTANTS AND RELATIONS 


In order to calculate a number of terms in the equations for 
the heat and kinetic energy transfers, it is necessary to know the values 
of a number of constants and to compute parameters that are functions of 
known variables. All temperature dependencies in the following rela- 
tions are written in terms of a temperature, T, measured in degrees 
Kelvin. 

1. The constants that are given below are found in the Smithsonian 


Meteorological Tables: 


R, = 287.05 J kg > oK + 
Babe UG1S5 1 oskgero Ree 
A -3 
Orie 1000. kg m 
OW] cleesn-boltmmenn constant) = 5.6687e100° J m- si> °KI 
a et ae | Be =) ook 
Serbian G7 To"Sany 4 ; or a 
2g Ser as gia 890) akg) (2k 


2. Since the tabulated emissivities of a water droplet are not readily 
available, the emissivity is approximated using the emissivity of a 
plane water surface. The emissivity of water, E, displays a weak 
linear dependence on water temperature. The linear dependence is 
ignored in this study for the range of droplet temperatures consid- 


ered and a value of E = 0.96 is used. 
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The approximation of the emissivity of a droplet by that of a plane 
water surface is not justified-when the wavelength of a radiation 
compares to a droplet radius. 

The dynamic viscosity of air as a function of temperature is given 


by (List, 1966, p. 394) 


2 B72 
u = 1.8325 x 1072} 226-16 + eo | iT | 


T + 120 


The units of viscosity are in Ky eae ia. 


The diffusivity of water vapor in air for temperatures between -40°C 


and 40°C is given by (Pruppacher and Klett, 1978, p. 413): 


e225) Cy | 1.94 


-4 
A = 0.211 x 10 | P [273.15 


The units of diffusivity are noest th, 


The specific heat capacity of water as a function of temperature is 


(Ibid, p. 89) 


coe 4186.34 [0.9973 Te Dee fone : 308,15)* 


ee dean ea bee 308.15)" forubiows (30K 


and, 


c = 4186.84 [1.0074 PPBe29) 8 9H0 (= 273.15)°| 
i form < 273; K 


-1 
The units for Cy are J se oK 


The saturation vapor pressure as a function of temperature is calcula- 


ted from the following polynomial in nested form (Lowe, 1977), 
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A, = 6984.505294, A, = 188.9039310 

A, = 2.133357675, oe ~1,2885809730x1072 
a 4.393587233x10 > aa ~8023923082x107° 
poe 6.136820929x107 14 


aan ite AgtT (Ay +T(Ay+T (A3+T (A, +T (A, tAGT))))) 


est 18 expressed in mb. 


The surface tension of water in air is (Pruppacher and Klett, 1978, 


p. 104) 


se = 10° 70k Oe SVEPUG: js 273.15)| FORT 1233" K 


The units of surface tension are J 7. 


The specific latent heat of vaporization as a function of temperature 


is expressed by the following empirical relation (Ibid, p. 89): 


h 
O78 ws 0) 


2, = 4186.84 [597.3 [223s (2) 
Vv 


pees 
: : -1 
where a is expressed in J kg ©. 
Alternatively, the specific latent heat of vaporization obtained 


by integrating Kirchhoff's Law is, 
2 = 2.501 x 10° - 2340.8(T ~ 273215) (b) 


The experimental relation for specific latent heat (Equation a) 
is used in the droplet thermodynamic and growth equations. Equation b, 
on the other hand, is used in the entropy flux equation of the ensemble 
model. The reason for using Equation b rather than Equation a is 


that the assumptions used in its derivation are then consistent with 
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some of the assumptions needed to derive the global entropy equation. 
A relation for the thermal conductivity of moist air is given by 
Beard and Pruppacher (1971) and by Pruppacher and Klett (1978, 

p. 418). The author wishes to indicate that the expression given 
by either reference contains typographical errors. 

The thermal conductivity of moist air can be given in terms of 
the thermal conductivity of dry air, Kae and of water vapor, ke 
The relation is obtained from the Mason-Saxena formula (see Bird et 
al, 1960, p. 258) for a gas mixture consisting of dry air and water 
vapor, where the amount of each gas is expressed in terms of its 
mole fraction (i.e. the ratio of the partial pressure of the 
gas to the total pressure). Since the mole fraction of water vapor 
in moist air is small, the thermal conductivity of moist air may be 
approximated by a first-order expansion of the Mason-Saxena formula. 


The result of the approximations is an expression given in terms of 


Ky and ko 


r 
k =k jt (1.17 - 1.02k /k4) | 


The constants in this relation are approximated using data for the 
viscosity of dry air and of water vapor. 
The temperature dependencies of Kg and kK, (Pruppacher and 


Klett, 1978, p.418) are given by, 


4.1868 x 10> 5.69 oli 273.15)] 
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The units of conductivity are given inJm = sS °K 
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APPENDIX H 


PROGRAM LISTING 


SIFSFTLISPPFRSEPSESSSPPSPEPRGSFPERPE LP RAPRSPEPPPSRPPREPSSPRKKSERSIDBDR 


THE DROPLET ENSEMBLE PROGRAM CALCULATES THE FLOW 
PROPERTIES OF AN ENSEMBLE OF DROPLETS THAT ARE 
ACCELERATING DUE TO THE AERODYNAMIC DRAG OF AN 
ACCELERATING MOIST AIR FLOW. UNITS USED ARE MKS. 
RESPSLELHELRSESLSESEPHESRLK ERATE PERE RHR SSE SHELSRKER RBH 
THE VARIABLES USED IN THE PROGRAM ARE GIVEN BELOW. 
NOTE THAT VARIABLES ENDING IN SUBSCRIPT 2 REFER TO 
CURRENT (UNKNOWN) VALUES. THE SUBSCRIPT 1 IS USED 

TO DESCRIBE VALUES OBTAINED FROM THE PREVIOUS 
CALCULATION. VARIABLES ENDING IN @ DENOTE VALUES 


tBe@ 
BED 
BEB 
BEB 


BEB 
£22 
REE 
EB 
BEB 
ES 


EE 
EES 
ES 
BES 
RES 
BRE 
EE 
EB 
EEE 
EES 
ERS 
Rs 
ERS 
EER 
BRE 
BRE 
£RE 
REE 
$s 
ERE 
ERE 
EES 
EEE 
SE 
ES 
eS *S 
ER 


SRS 
SE 
EBB 
BES 
EES 


EES 
ERE 
EBS 
ERE 
ESE 
HEB 
EER 
£82 
BRS 


THE PLENUM. 
VARIABLES 
X AXIAL DISTANCE FROM PLENUM 
AREA FLOW CROSS-SECTIONAL AREA 
VELX ONE-DIMENSIONAL VELOCITY OF MOIST AIR 
TD TEMPERATURE OF MOIST (AND DRY) AIR 
P PRESSURE OF MOIST AIR 
UC ONE-DIMENSIONAL VELOCITY OF DROPLETS 
TC TEMPERATURE OF DROPLETS 
TVIRT VIRTUAL TEMPERATURE OF MOIST AIR 
RAD DROPLET RADIUS 
ROD DENSITY OF DRY AIR 
RHO DENSITY OF MOIST AIR 
N NUMBER OF DROPLETS PER UNIT VOLUME 
ROS MASS OF WATER DROPLETS PER UNIT VOLUME 
RMIX MIXING RATIO OF MOIST AIR 
SPH SPECIFIC’ HUMIDITY OF MOIST AIR 
RH RELATIVE HUMIDITY OF MOIST AIR 
FD MASS FLUX OF DRY AIR 
FA MASS FLUX OF MOIST AIR 
FV MASS FLUX OF VAPOR 
FS MASS FLUX OF SPRAY (DROPLETS) 
FW MASS FLUX OF WATER SUBSTANCE 
ES SATURATION VAPOR PRESSURE OF MOIST AIR 
EA VAPOR PRESSURE OF MOIST AIR 
FNUM DROPLET NUMBER FLUX 
MASS DROPLET MASS 


THE FOLLOWING CONSTANTS, DERIVED PARAMETERS AND 
VARIBLES ARE REQUIRED IN THE CALCULATIONS : 


CONSTANTS 


STEFAN-BOLTZMANN CONSTANT 
EMISSIVITY OF WATER 
DENSITY OF WATER 


:GAS CONSTANTS FOR DRY AIR AND 

:WATER VAPOR, RESPECTIVELY 

:RESPECTIVE SPECIFIC HEAT CAPACITIES AT 
:CONSTANT PRESSURE 

:RESPECTIVE SPECIFIC HEAT CAPACITIES AT 
:CONSTANT VOLUME 
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pte DERIVED PARAMETERS AND VARIABLES 

riecieg VIS VISCOSITY OF AIR 

bed RE REYNOLDS NUMBER 

liga NU NUSSELT NUMBER 

hts SH SHERWOOD NUMBER 

fi PR PRANDTL NUMBER 

gett SC SCHMIDT NUMBER 

ndtalig DIFF DIFFUSIVITY OF MOIST AIR 

edi TENS SURFACE TENSION OF WATER IN AIR 

dg CURV FACTOR INCREASING THE DROPLET SATURATION 
bc ed VAPOR PRESSURE DUE TO CURVATURE EFFECT 
‘ghd K THERMAL CONDUCTIVITY OF MOIST AIR 

th dg CP SPECIFIC HEAT CAPACITY AT CONSTANT 

finda) PRESSURE FOR MOIST AIR 

72% CV SPECIFIC HEAT CAPACITY AT CONSTANT 

ified VOLUME FOR MOIST AIR 

bendid VHEAT SPECIFIC LATENT HEAT (EMPIRICAL 

Faded RELATION USED IN DROPLET EQUATION) 


SREESSTERPERPESSRPSELPFERKREFERL HSER ETEK LSE KRKRLPELSKEFSEFRSKRESSVSS 


##**%®® CALCULATIONS FOR SPRAY AND AIR FLOW #######48828884 
##*88*® PROPERTIES BEGINNING AT THE SPRAYERS. #*#######3244 
#** FIRST, CALCULATE THE MOIST AIR PROPERTIES 
*** AT THE SPRAYERS. THE CALCULATIONS ASSUME FRICTIONLESS 
#**® ADIABATIC FLOW (DETERMINED BY INITIAL CONDITIONS 
*** AT THE PLENUM). PART TWO OF THE PROGRAM 
*** CALCULATES THE COUPLED FLOW OF AIR AND SPRAY. 
#** A VARIABLE GRID SPACING IS USED. 
#** X1 IS THE SPRAYER LOCATION. XMAX INDICATES THE 
##* BEGINNING OF THE WORKING SECTION. 
PERPRKRSERSKFKRPERPKERP RRS KKK EK SREREKEKLLERRPREKRLEKRKERERPRRRKRRD 
IMPLICIT REAL*8(A-H,0-Z) 
REAL SCALE*8(3) ,CFIG1 (10) ,CFIG2(10) , TITLE(19) 
REAL*8 PI/3. 1415926535898/ ,SIG/S.6687D-08/, 
1 EMIS/@.96D0/,DENS /1000.DO/ 
REAL*8 RD/287.05D0/ ,RV/461.51D@/,CPD/1005.D0/, 
1 CPV/1850.D0/, CVD/718.DO/,CVV/1390.DO/ 
REAL*8 X2, AREA2, RMIX2,P2,RHO2, VELX2,TD2,UC2,TC2 
REAL#8 N,MASS,NU, LATENT ,K 
READ (5,10) P@,TDO,UO, RH 
10 FORMAT (E12.5,F7.2,F8.3,F7.2) 
READ(5,20) RAD,TC1,X1,UC1,ROS 
20 FORMAT (E10.3,F8.3,F7.4,F8.3,F7.2) 
KERPKESEEKESEKRREFEEREP EERE RKP KEE EEK ERRRERBPELSRRSKRRSKEPESRD 
*#* THE INCREMENTATION SCHEME DESCRIBED BELOW 
##* APPLIES FOR THE GEOMETRY OF THE NRC TUNNEL. 


*%*® XMAX = BEGINNING OF THE CONTRACTION SECTION 


##* THE CONTRACTION DISTANCE IS DIVIDED INTO 
##® 3 SECTIONS. THE FIRST TWO SECTIONS USE 

##* LOGARITHMIC INCREMENTS.THE INCREMENTS ARE 

##® DETERMINED FROM THE RELATION, 

#8 DELTA = (LOG(XTWO) - LOG(XONE))/N 

##® WHERE XONE, XTWO ARE ENDPOINTS OF THE SECTION 

#** CONSIDERED AND N IS THE NUMBER OF INCREMENTS 

**® IN THE SECTION. THE THIRD SECTION USES CONSTANT 

##* INCREMENTS DETERMINED BY THE RELATION, 

eae DELTA = (XTWO - XONE)/N 

##* DELI, DEL2, DEL3 DENOTE THESE INCREMENTS. 

ke 

#*® XC1, XC2 IN THE PROGRAM CORRESPOND TO XONE, 

##® XTWO OF THE FIRST AND SECOND SEGMENTS. 

##* XOUT SHIFTS OUTPUT FROM EVERY .05 M. TO .10M. 
LERDREREREESREEELELSERESER ELE LSE SELELRELREKEE SERS RLVS BR 
#%* THE MODEL PARAMETER °EPS’ DENOTES THE ALLOWABLE 

##® RELATIVE ERROR AND °’LIM’ DENOTES THE MAXIMUM NUMBER 
##% OF ITERATIONS ALLOWED IN SOLVING FOR VELOCITY. 
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READ (5,30) XMAX,N1,N2, N3 LIM, EPS, EPS2 
30 FORMAT (F7.4,315,14,2E8. 

READ(5,40) XC1,XC2, XCHG. 

40 FORMAT (3F9.6) 

RADOUT=RAD/ | .@E-06 
SEPRETSPRESSPSPKE SPER PERLE RSLTSLERTEPRRSRLELELTERRFETPRRSPEKERRPRSKE 
*#*® PROGRAM PRINTS DESCRIPTION OF TUNNEL. ARRAY 
##® °TITLE’ GIVES TUNNEL IDENTIFICATION. ARRAYS 
*#® °CONFIG1’, CONFIG2’ DESCRIBE TUNNEL GEOMETRY. 


SEES PLE SERS RETR SKKEPEKE SSS REE KRKKLHKE KKK KKK ESET KKK RS KRKE RK BRB 


READ (3,50) (TITLE(I), I=1, 10) 
50 FORMAT (10A4) 
READ(3,60) (CFIGI(1),I=1, 10), (CFIG2(I),I=1, 10), 
1 (SCALE(1) , 1=1,3) 
60 FORMAT (10A4/10A4/3(A8/) ) 
DEL 1 = (DLOG(XC1/X1))/N1 
DEL2= (DLOG (XC2/XC 1) ) /N2 
DEL3= (XMAX-XC2) /N3 
MASS=4./3.®P1*RAD*RAD*RAD* DENS 
N=1.@D-03*ROS/MASS 
AREAO=CS (@. QODOO) 
E=RH*VAPOR (TD@) 
PD=P0-E 
RMIX1=0.62198*E/PD 
SPH=RMIX1/(RMIX1+1.) 
TVIRT=TDO* (1. +SPH#0.6077688) 
CP=CPD* (1. +(CPV/CPD- 1.) *SPH) 
CV=CVD* (1. + (CVV/CVD-1.) *SPHD 
GAM=CP/CV 
ROO=PO/(RD*TVIRT) 
FAQ=ROO* AREAQ*UO 
B1=2.0*GAM*PO/ (ROO* (GAM-1.)) 
B=B1# (UQ* AREAQ) *# (GAM-1.) 
WRITE(6,70) TITLE, XMAX,CFIG1,CFIG2,TDO,P,UO,RH 
7@ FORMAT (?1’,T12, THERMODYNAMIC PROPERTIES OF SPRAYS’/ 
> + \T13,’AND AIR FLOW IN AN ICING TUNNEL’ / 
"9 TS, TUNNEL DESCRIPTION? / 
, ding SS SSeS Soe sien 
°@’ ,T9, TITLE: ’/ 
9 T12, 10A4/ 
> + °T9, CONTRACTION LENGTH: */ 
*°5°T12,FS.4,’ METERS’ / 
> > TQ, CONFIGURATION: ?/ 
"eS" TI21OA47*  * 12, 1OA4/ 
°@’ T9, TUNNEL OPERATING CONDITIONS AT BELL MOUTH’ 
/ > *\T12,’ TEMPERATURE: ’ ,F8.2, 
SS DECEC YY Maia yT iz: 
*PRESSURE:’ ,F9.0,’ PASCALS’/ 
* + T12,’ AXIAL VELOCITY:',F9.3,’ M/,SEC’/ 
* 9°T12)RELATIVE HUMIDITY’ ,FS.2,’ 2°) 
AREA1=CS(X1) 
U=AREAO*UO/AREA 1 
PRELKEEHKKKEKKKEFKERKERKRPE LEEK EKERPESP FREE FEEL RLERBKEREREEKEPRS 
#***% LOOP TO SOLVE FOR MOIST AIR VELOCITY *#########29 
#** IF RELATIVE CHANGE IN VELOCITY ESTIMATE IS 
#** GREATER THAN THE ALLOWABLE ERROR EPS, CONTINUE 
#** ITERATIONS. IF NUMBER IS GREATER THAN LIM, PRINT 
#*® ESTIMATE, X2, RESIDUAL, NUMBER OF ITERATIONS 
#*® AND TERMINATE THE CALCULATIONS. 
KRERREEERKKKEKKEKREKEKEER KSEE PEKSKFKKSEPKKKPRSSERSERRREB 
DO s@ J=1,LIM 
U1=((UO*UQ-U*U+B 1) / (B/AREA 1## (GAM-1.)))*#(1./(1.- 
1 GAM)) 
ERR 1=DABS( (U1-U) /U1) 
FN=U 1#U1+B* (AREA1#U 1) ## (1.-GAM) -U0*U0-B1 
KLIM=J 
IF (ERR1.LE.EPS) GO TO 90 
U=U1 
8@ CONTINUE 
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GO TO 340 
PRERTERLEPKEKEREPRETRKKRE PE PSR RRKKLHKERRFERRRPRSERPRTETEBRLSTPEBDR 
##® SOLVE FOR REMAINING AIR AND DROPLET 
#* PROPERTIES AT SPRAYERS. 
SESPPRSSSPEPERPEPE SSSR RESETS ISDE SE SSPE ETFLERPERSREPFERKRESERSSSRSP 
98 CONTINUE 

VELX 1=U1 

RHO 1=FA@/ (AREA 1*®VELX 1) 

ROD=RHO1/(RMIX1+1.) 

P 1=PQ* (RHO1/ROQ) **GAM 

TVIRT=P 1/(RD*RHO1) 

TD1=TVIRT/(1.+0.6077688*SPH) 

VIS=1.8325D-05* (296. 16+ 120.)/(TD1+120.)® (TD1/296. 16) #® 


Lies 
SSRELFESESEPRRPPESKKEPERPSRKPPKEREPTPRELLETFSRISERESSREREBPSE 


#?% IF INJECTION VELOCITY IS CODED -1.@, DROPLETS 

*#® MOVE WITH THE AIR VELOCITY. IF INJECTED 

#*#2® TEMPERATURE IS CODED 0.0, DROPLETS ARE INJECTED 

*** WITH AIR TEMPERATURE. 

KESPEKRERSL ARERR PK SSA LTAEREPRESPPEKKS RPE SREPSERRKE RSLS BRR 
IF(UC1.EQ.-1.0) UCI=VELX1 
IF(TC1.EQ.0.00) TC1=TD1 


EKER RE PKR EPKETE REE PERK KEPRREFE PEPER BR PRKSKL ERE TETKRELRRERKSE 


***® CALCULATE DRAG COEFFICIENT 
KPEFRLKKERKK KERR ERE KKK REKKRSREKKSERPKREPRKERESRRSKPEREKPSEELSR 
RE=2. *RAD#*RHO1* (VELX1-UC1)/VIS 
ALPH=0. 189D00 
BETA=0.632D00 
IF(RE.LT.20) GO TO 100 
GO TO 120 
160 CONTINUE 
IF(RE.GT.1.5) GO TO 110 
ALPH=@. 102D00 
BETA=0. 955D00 
GO TO 120 
110 CONTINUE 
ALPH=0. 115D00 
BETA=0. S02D00 
126 CONTINUE 
C1=DENS*CW(TC1)#RAD*RAD/3. DOO 
DRAGF= 1. +ALPH®RE#*BETA 


HFKERKEREREPESERE SEEK KEKTE KEPT F KR KKETRFESRPRPRESEPRPPREPKRSSR PRS 


#*® ACCOUNT FOR VENTILATION EFFECTS 
REKKEKKPPRK SRST SRPDREKKKE PRESET RPK SE RERKFEPERPSRERPSEVRKEREKRS 
TEMPF=(TD1+TC1)/2. 
DIFF=0.211D-04* (101325. /P1)# (TEMPF/273. 15) #*1.94 
PR=CP#VIS/COND(TD1 , RMIX1) 
SC=VIS/(RHO1*DIFF) 
PR1=PR**(1./3.) 
SC1=SC*## (1./3.) 
RESRT=DSOQRT (RE) 
XVAR=SC 1#RESRT 
YVAR=PR1#RESRT 
IF (XVAR.GE.1.4) SH=2.*(@.78+0.308*XVAR) 
IF(XVAR.LT.1.4) SH=2.*(1.00+0. 108®XVAR*XVAR) 
IF(YVAR.GE.1.4) NU=2.*(0.78+0.308*YVAR) 
IF(YVAR.LT.1.4) NU=2.*(1.00+0. 108*YVAR*YVAR) 
EA=P 1*RMIX1/(@.62198+RMIX1) 
TENS= 1 .OD-03* (76. 10-0. 155* (TC1-273. 15)) 
CURV=DEXP (2. *TENS/ (DENS*RV#TC 1*RAD) ) 
ES=CURV# 100. *VAPOR(TC1) 
EF= (ES+EA) /2.D00 
VHEAT=LATENT (TC1) 
F 1=VHEAT*DIFF*SH* (ES-EA) /(2. *RV®TEMPF* 
1(1.-EF/P1)) 
K=COND(TD1,RMIX1) 
F2=K*NU* (TC1-TD1)/2.DO0 
TRAD=TC1#TC1*TC1*TC1-TD1*TD1*TD1*TD1 
F3=SIG#EMIS*RAD*TRAD 
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BE 1S2.*PI*RAD*DIFF® (ES-EA) *SH* VHEAT/ (RV*TEMPF® (1 . -EF/ 
FF2=2.*P1*RAD*K*NU* (TC1-TD1) 
FF3=4.*P1# (RAD*RAD) *SIG*EMIS*TRAD 
DQ=F 1+F2+F3 
DTD=DQ/C1 
SSSELERPEETELEPESRSEREK SPE KTLKREPRPKTERKRERREKRRKFFLKEREKESBB 
##® CALCULATE TERMS OF THE GLOBAL ENTROPY EQUATION AT 
**#® THE SPRAYER POSITION. ENTROPY TERMS ARE EHUM, EPRESS, 
##® ETC. NOTE THAT THE SPECIFIC LATENT HEAT IS 
#*® CALCULATED FROM KIRCHHOFF’S LAW. 
EPEEEL GREE HEL EEE TERE KE SE RL ESLER KRERRE PRESS SRORE ERED 
ROS= 1 .@D-03*ROS 
FS=ROS*UC1*AREA | 
RH=EA/(100.*VAPOR(TD1)) 
RHPER=RH* 100. 
FS0=ROS* AREA 1#UC1 
FD@=ROD* AREA 1*VELX 1 
FV=RMIX1#FDO 
FWO=FSO+FV 
FNUMO=N*UC1*AREA1 
CPW=4218.D00 
COEFF=FV*CPW+FD0*CPD 
IF(RH.GT.1.@D-03) GO TO 130 
EHUM=0. 00D 
GO TO 140 
130 CONTINUE 
EHUM=FV*RV*DLOG (RH) 
140 CONTINUE 
EPRESS=FDO* RD*DLOG (P 1/ (RMIX1/0.62198+1.)) 
ELAT=FV#*SPLAT (TD1)/TD1 
ESENSC=FS*CPW*DLOG (TC 1) 
ETEMP=COEFF*DLOG (TD 1) 
WRITE(6, 150) ROS,RADOUT,TC1,TD1,UC1,VELX1,RHPER 
15Q FORMAT(’@’,T9,’SPRAY PROPERTIES AT NOZZLES: °/ 


1 etl iy LOUID: WATER, CONTENT :2 (FG. 2, 

ie  PKGOM/ CUCM 7 2 1125 MEAN DROPLET (SIZE?.? 

3 seO.O 70 UM ICRONS 7 9) 2, SPRAY) TEMPERATURE: 
4 SES 525% DEG C7 ak, slile:: "AIR TEMPERATURE: ’ 

2 whode)  URG=C yan ba ee SORRAY VELOCITY :” 

6 AO co DAMA Stee nn Tio “AIR VELOCITY: 

7 ME icce EM SOG Is | an eT i2. 

8 -RELATIVE HUMIDITY: ’ ,F8. ie Nee 


WRITE(G,1600) XT SNF, SCALE (1), DELI, XC1 NZ, SCALE(2) ,DELZ2, 

1XC2,N3,SCALE(3) , DEL3 , XMAX 
160 FORMAT (°@” , T9, >X-RANGE? ,T22, ’ NUMBER’ , T31, *SCALE’ ,T41, 
“INCREMENT #7 oe 1 F8.6, -/7? a PKG): att 122) 
I5,129,As,1T4@, E13. sae OP es SG Ses 6, 

oh 


Sabie ting 122, PS 1295 Ae, 40,n 10.77, 
6, vy, 


8. 
ET iG, re 122, [5,129 AS 40, lo. 7, 
: TO" F8.6) 
WRITE(6, 170) 
170 FORMAT @ 0? 5 T4,°X", TOS" AIR: SREED” | TZ21 5° DROP SPEED*.,TS1, 
MAUR TTEMP], fiat. UROF, TEMP” ,TS35 PRESS’ ,T61, 
UDENS PLY 3170, BREA o, 17S, “ERROR 191, 
PPT ERe/4 
Se 1S Ro BRA 4 bc, UWC > TS6," RADIUS”, 
T48, ’NUMBER DENSITY’ 7//) 
WRITE(6, 320) X1LSVELX1,UC1,TDL,TC!,P1,RHO1, AREAT,ERR1, 
IKLIM, RHPER, RMIX1, ROS, RAD, N 
FENTO=ESENSC+ETEMP- -EPRESS+ELAT-EHUM 
TD2=TD1 
XOUT=9.2SD0 


ESREFEFPREFSEREKSSRHSPLEPFPSRESPTESRREKLERPSLPRSRSTPRPPPRPPRRETD 


*#*® PART TWO OF PROGRAM CALCULATES PROPERTIES OF 
*%® COUPLED AIR AND SPRAY FLOWS. 
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198 
200 


DO 330 ICOUNT=2, 16000 


IF(X1.GT.XC2) GO TO 18@ 
IF(X1.LT.XC2) DEL=DEL2 
IF(X1.LT.XC1) DEL=DEL1 
X2=DEXP (DLOG(X1) +DEL) 
DINCR=X2-X1 

GO TO 198 

CONT INUE 

X2=X1+DEL3 

DINCR=X2-X1 
IF(X2.LT.XMAX) GO TO 190 
X2=XMAX 

DINCR=XMAX-X1 

XOUT=X2 

CONTINUE 

AREA2=CS (X2) 

CONT INUE 


C SFPSRFREPPELRERE SKK PSRKRPKPRERKLPRBRERSRSESPRERTRSLGRERSPLRLSSD 


C **#® CALCULATE SPRAY PROPERTIES 


C FFERPSEPR KES RRPEKLEKKEKSEK KEK EPL RKEPSSTERPRTRIRTSRRRSRBSSS 


210 
220 


Q=(9.*VIS) / (RAD* RAD* DENS) * DRAGF 
UC2=DSQRT (UC 1#UC 1 +Q*DINCR# (VELX1-UC1)) 
IF(UC1.LE.1.@E-06) GO TO 210 

RAD= (RAD* RAD-2. *F 1*DINCR/ (DENS*# VHEAT#UC1) ) ##0.5 
TC2=TC1-DTD*DINCR/UC1 

IF (X2.EQ.XMAX) GO TO 310 

GO TO 220 

CONT INUE 

TC2=TC1 

CONT INUE 

MASS=4.73.*P 1*RAD*RAD* RAD* DENS 

N=FNUMO/ (AREA2*#UC2) 

FNUM=N*UC2# AREA2 

ROS=MASS*N 

CMOM 1=FS*UC1 

FS=ROS* AREA2*UC2 

CMOM2=FS*UC2 

CMOM=CMOM 1-CMOM2 

U=VELX 1* AREA 1 /AREA2 

RMIX2= (FWO-FS) /FDO 
COEFF=FV*CPW+FDO*CPD 
TENSsb..OD>5037 ©76.510-0. ISS? (TCZ2-273..15)) 
CURV=DEXP (2. *#TENS/ (DENS*RV*TC2*RAD) ) 
ES=CURV* 100. * VAPOR (TC2) 

FV=RMIX2*FDO 

SPH=RMIX2/ (RMIX2+1.) 
PCON=RMIX2/0.62198+1. 
TCON=1.+0.6077688*SPH 

AVCS= (AREA1+AREA2) /2. DO 

FORCE=6. *N*PI*VIS*RAD*? (VELX1-UC1) *DRAGF*® AREA 1*DINCR 


C FPFPFRSPPRPSK RFE PK SKE SSSRARERREKEKPREPRPRSRRER KSEE RPSL SSRPSD 


C *#** LOOP TO SOLVE GLOBAL EQUATIONS FOR AIR PROPERTIES 


C FPSPPRPERHS EPH PRSSFE SFR REFEREE ELPSEREERSKERSRRSRKERKE LES SSS 


230 
Cc 


DO 230 KK=1,LIM 
UHOLD=U 
PTERM=FD0* ( (VELX 1 -UHOLD) + (RMIX1*VELX 1-RMIX2* UHOLD) 
) +CMOM 
P2=P 1+ (PTERM-FORCE) /AVCS 
EPRESS=FD0O* RD* DLOG (P2/PCON) 
ELAT=FV*SPLAT (TD2) /TD2 
ESENSC=FS*CPW* DLOG (TC2) 
FENT =-ESENSC+EPRESS-ELAT + EHUM 
YT2= (FENT +FENTO) “COEFF 
TD2=DEXP (YT2) 
TVIRT=TD2*TCON 
U=FDO*RD*#TVIRT* (1.+RMIX2) /(P2*AREA2) 
ERR2=DABS ( (U-UHOLD) /U) 
KLIM=KK 
IF (ERR2.LE.EPS2) GO TO 240 

CONTINUE 


GO TO 360 
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#*® CALCULATE QUANTITIES CHARACTERIZING DROPLET 
**® FLOW RATES, PARAMETERS, ETC. THAT ARE FUNCTIONS OF 
**® AIR PROPERTIES. THESE RESULTS ARE USED IN THE 
##® DROPLET FINITE-DIFFERENCE EQUATIONS TO SOLVE FOR 
#*#® UNKNOWN PROPERTIES AT NEXT INCREMENT. 
SREDERTEKETRL REPRE KKK PEER LPKTKKE SK KEKEKKERERPKERRSPERPRHRRBREM 
240 CONTINUE 
VELX2=U 
EA=P2*RMIX2/(@.62198+RMIX2) 
RH=EA/ (100. *VAPOR(TD2) ) 
IF(RH.GT.1.0D-03) GO TO 250 
EHUM=9 . ODO 
GO TO 260 
250 CONTINUE 
EHUM=FV*RV*DLOG (RH) 
268 CONTINUE 
VIS= 1 .8325D-05* (296. 16+ 120.) /(TD2+ 120.) * (TD2/296. 16) 


RHO2=FDO* (1. +RMIX2) / (AREA2*VELX2) 
RE=2. *RAD*RHO2* (VELX2-UC2)/VIS 
ALPH=@. 1S90000D00 
BETA=0.632D00 
TECRES LT .26). GO" TOe 270 
GO TO 290 
270 CONTINUE 
IF(RE.GT.1.5) GO TO 280 
ALPH=0. 102D00 
BETA=0.95SSD00 
GO TO 290 
280 CONTINUE 
ALPH=0. 115D00 
BETA=0.802D00 
290 CONTINUE 
C1=DENS*#CW(TC2) *RAD*RAD/3. DOO 
DRAGF= 1 .+ALPH*RE**BETA 
300 CONTINUE 
PR=CP*VIS/COND (TD2, RMIX2) 
TEMPF= (TD2+TC2)/2. 
DIFF=0.21 1D-04* (161325. 7P2)* (TEMPF/273. 15) °* 1.94 
SC=VIS/ (RHO2* DIFF) 
FREER ON 732 
SC 1=SC#*(1./73.) 
RESRT=DSQRT (RE) 
XVAR=SC 1*RESRT 
YVAR=PR1*RESRT 
IF (XVAR.GE.1.4) SH=2.*(@.78+0.308*XVAR) 
IF(XVAR.LT.1.4) SH=2.*(1.00+0. 108*XVAR*XVAR) 
IF(YVAR.GE.1.4) NU=2.*(0.78+0.308*YVAR) 
IF(YVAR.LT.1.4) NU=2.*(1.00+0. 108*YVAR*YVAR) 
EF=(ES+EA) /2.D00 
VHEAT=LATENT (TC2) 
F 1=VHEAT*DIFF*SH? (ES-EA) / (2. ®RV*TEMPF® 
Ll bezbe)) 
K=COND (TD2 , RMIX2) 
F2=K*NU* (TC2-TD2) /2. DOO 
MRAD=fe2" 1C2*TC2* Te2—-1D2STD2°Tp2°T bz 
F3=SIG*EMIS*RAD*TRAD 
DQ=F1+F2+F3 
DTD=DQ/C1 
TD1=TD2 
TC1=TC2 
RMIX1=RMIX2 
UC 1=UC2 
VELX 1=VELX2 
X1=X2 
P1=P2 
RHO1=RHO2 - 
AREA 1=AREA2 
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IF (X2.LT.XOUT) GO TO 330 
XOUT=XOUT+0.0S 
IF (XOUT.GT.XCHG) XOUT=XOUT+0.95 
3190. . CONTINUE 
FF 1=2.*P1*RAD*DIFF* (ES-EA) *SH* VHEAT/ (RV*TEMPF®? (1.- 
lL ERZP2)) 
FF2=2.*P1*RAD®K*NU* (TC1-TD1) 
FF3=4.*P1* (RAD*RAD) *SIG#EMIS#TRAD 
RHPER=RH* 100. 
WRITE(6,320) X2,VELX2,UC2,TD2,TC2,P2,RHOZ,AREA2, 
I GERR2Z, KLIM, RHPER , RMIX2,ROS,RAD,N 
320 FORMAT(’ Q°, F6.4,4F 10. 4, F507 2F8.4,E13.4, 16/ 
1 : > F5.1,E13.4,F11.6,2E13.5) 
IF (X2.EQ.XMAX) STOP 
330 CONTINUE 


STOP 
PEDELERLTEPEPRTEPPELLELPKLESEPEPKRERERPKKE SPER SPRSPEPRERPKLIPSTESE 


**® PRINTOUT IS GIVEN WHEN ITERATIONS REQUIRED TO 
##® SOLVE VELOCITY EQUATION ARE GREATER THAN LIM. 
SRSPETFSPLKREKPPESE KEKE KEK RPK SKEKERSEKKEDKE SEK EKTERELRPRERKSRSS 
348 CONTINUE 
WRITE(6,350) U1,X2,FN,J 
350 FORMAT (’-? ,TS,'U= *, F8.4,* DOES NOT YIELD ROOT OF EQUATION? 
1 / ,TS,’X= °,F10.6,T17,FUNCTION= °,E10.2, 
2 T40,’NO. OF ITERATIONS =’, 14) 
ey Olas al WUE ce ek ea 
#9 PRINTOUT IS GIVEN WHEN ITERATIONS REQUIRED TO SOLVE 
##® GLOBAL EQUATIONS ARE GREATER THAN LIM. 
RRRPPREEFKEKRESEPTE FERRER KEKKREEKK KT EKEPESPE KP SEK PPRKKEBRS 
360 CONTINUE 
WRITE (6,370) _VELX2,X2,KLIM 
370 FORMAT (?-*,TS,’U =',F12.7, DOES NOT SOLVE SYSTEM’/’ °, 
1 T550X Ls PIO! 6,T23,’NO. OF ITERATIONS =? 
2 14) 
STOP 
END 


PEBEKREERLKEKREPRERPEKFEERREREERESRRKEFEKREEPRRERRKRERLKLEKERRR RHE 


#*® SUBPROGRAM TO CALCULATE TUNNEL CROSS-SECTIONAL AREA. 
KERREEEKREF ERR SRKEREKP KEEP REEKKERRERKRREKEKKEKRKEKEEKER KEKE 
REAL FUNCTION CS*8(X) 
REAL*#8 X,XHOLD1! , XHOLD2 
XHOLD 1=X-@. 1524 
IF(X.GT.0.1524) GO TO 10 
XHOLD2=0 . 38 10- (@ .02323-XHOLD1*XHOLD1)##0.5 
CS=4.*XHOLD2*XHOLD2 
RETURN 
10 CONTINUE 
XHOLD2=0 . 2286- (0.03846*XHOLD1) 
CS=4 .. *XHOLD2* XHOLD2 
RETURN 
END 


RKBKERERKKKSREKPREKERKKEK RE RP KE RPESEKRKELEPRPR RETR RRR RRP 


*#% SUBPROGRAM TO CALCULATE SATURATED VAPOR PRESSURE. 
ESSERE SESE SESE ESE SKELTER SEAR RK KEE SRE REE R ERE 

REAL FUNCTION VAPOR*8(T) 

REAL*8 AQ,A1,A2,A3,A4,A5,A6,T 

AQ=6984 .SO0S5294 

Al=-188.9039310 

AZ=2 133357675 

A3=-1.2885809730D-92 

A4=4.393587233D-95 

A5=-8 .023923082D-08 

A6=6. 136820929D- 11 

VAPOR=AO+T# (AL+T* (A2+T# (A3+T® (A4+T® (AS+A6*T) ) ))) 

RETURN 

END 
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*#® SUBPROGRAM TO CALCULATE LATENT HEAT OF EVAPORATION 
**® USING AN EMPIRICAL RELATION. RESULT IS USED IN 
##* CALCULATIONS FOR DROPLET EQUATIONS. 
PKEPESESSEERRE SPE EKKE RET KK KSKREEPEESPRL SER RSRSEKEBRREKRKRERP 
REAL FUNCTION LATENT*S(T) 
REAL*§ T 
EXP=0. 167+3.67D-04*T 
LATENT =597 .3* (273. 15/T) #*EXP 
LATENT=LATENT#4 186.84 
RETURN 
END 
PSRVEETEEEPERPKLEKFPEEKEPETKPE PERT FP PREP KPERRRRERPFPPRRPKPREKF 
##® SUBPROGRAM TO CALCULATE LATENT HEAT OF EVAPORATION. 
##® RELATION IS DERIVED FROM KIRCHHOFF’S LAW AND IS 
#*#® USED IN ENTROPY EQUATION. 
RERPVKEKKEKKREPFFRRPKERPEKEKKKVRERFP PKR RLPTESRETEREEERPEPERBRBD 
REAL FUNCTION SPLAT#8(T) 
REAL*S T 
SPLAT=-2340 . S@DOO* (T-273. 15) +2.501D06 
RETURN 


END 
PSELERSEEKKEPPRSPP TPT ERE SELERLSPRREPFPRPETRPRPEPRELPBVESPBSEBS 


##® SUBPROGRAM TO CALCULATE SPECIFIC HEAT CAPACITY 
#*#® OF WATER. 
SERRELETFFERPERPEREFSEKEKESERPFERPKREEFKLSREKEAPRPEEPPREPRELPRESBPPS 
REAL FUNCTION CW*8(T) 
REAL*8 T,T1 
IF(T.LE.273.15) GO TO 10 
T1=T-308. 15 


CW2059979+3 710-06? f1*11tS.SD-O9eT Terie hier t 
CW=CW*4186.84 
RETURN 

1@ CONTINUE 
CW=1 5007448. 29)-05? Cl =273.2'5) * 1-2 7/3215) 
CW=CW*4 186.84 
RETURN 


END 
PKBSKEKKLKKESKRPSEPSE SRK KP KEKE PKR KELSKRRSSRPSPEPPEERKVSETERRESD 


*#® SUBPROGRAM TO CALCULATE THERMAL CONDUCTIVITY 
##* OF MOIST AIR. 
REPRE KRKELERKRERECKERSRSKKE KK KEKKRKKETRKREEEKLHSRKVKEKRSRERSKEER 
REAL FUNCTION COND#8(T,RMIX) 
REAL*§ T,RMIX,EFAC 
DCOND=4. [S6D-03* (5.69+0.617* (T-273. 15)) 
VCOND=4. 186D-03* (3.78+0.020# (T-273. 15) ) 
EFAC=RMIX/ (@.62198+RMIX) 
COND=DCOND* (1.-(1. 17-(1.@2*VCOND/DCOND) ) *EFAC) 
RETURN 
END 
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